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PREFACE 


The main objective of this study is to develop a numerical model 
for simulating the movement of well graded sediment through a stream 
network. In Part 2 of this report the theoretical background of the 
model is described and the governing equations are formulated. The 
numerical solution of these equations is presented in Part 3., Hydraulic 
routing can be performed using any acceptable algorithm supplied by the 
user because the water movement is assumed to be uncoupled from the 
sediment process. The model can be used in conjunction with any 
suitable sediment yield model to supply the water and sediment runoff 
from lateral areas. The application and results of the model are 
discussed in Part 4. The computer program for the numerical model of 
the East Fork River system is described and listed in the third 
addendum. 
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_ To convert _ 
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INTRODUCTION 


'SI 

This report describes a one-dimensional numerical model designed to 
simulate sediment transport in natural channels. The physical processes 
associated with the sediment movement are reproduced using a variety of 
algorithms. These algorithms incorporate sets of equations that operate 
on input data in a predetermined sequence to generate output data repro¬ 
ducing the actual physical process. A satisfactory model must include 
the more relevant aspects of that process. 

Sediment moves driven by hydrodynamic forces exerted by the flow of 
water which in many instances is highly time dependent. The sediment 
transport model must, therefore, account for unsteadiness in sediment 
movement. 

The dependence of sediment motion on flow conditions makes it also 
dependent on the longitudinal variations the flow experiences as a 
result of stream boundary irregularities. These variations are 
reflected in the spatial variability of the sediment load distribution. 
Depending on particle size, some particles may be carried primarily in 
suspension, while others move entirely as bed load. In addition, 
depending on flow conditions, particles moving in suspension at one 
place may be moving as bed load farther downstream. ^ Whether a 
particular size fraction moves primarily as suspended load or bed load 
determines to what extent that fraction of the sediment load will lag 
behind the flood wave and therefore determines what the magnitude of 
longitudinal sorting will be. This means that the transport model must 
reflect the dependence of the sediment load lag on the material 
properties of the sediment as well as on hydraulic conditions. Whenever 
the bed material consists of a mixture, its transport involves the 
motion of a multitude of particles of diverse sizes. Some particles may 
deposit on the streambed while others are scoured away, resulting in a 
size composition of the material in transport different from that of the 
bed. A realistic model must account for the interchange between the bed 
surface material and the moving sediment load, and should simulate the 
residual transport capacity of the stream. The latter is a measure of 
the ability of the flow to further entrain material of a given size 
fraction in the presence of all the fractions already in motion. 


During the above particle interchange, the bed material particle 
size composition changes continuously and, in the process, the bed may 
experience a net amount of aggradation or degradation. For certain flow 
conditions the bed degradation in a reach may be limited by the 
formation of an armoring layer, over which sediment may move either in 
suspension or as intermittent bed-load waves. The armoring layer may be 
destroyed during high flows and reformed at subsequent low flows. A 
model must therefore be capable of tracking the streambed profile 
evolution and the changes in bed material size distribution. 

The proposed model is designed to meet the above criteria. It can 
simulate the unsteady transport of sediment mixtures through a network 
of nonbifurcating channel reaches, and it can be used in tandem with a 
suitable sediment yield model, like the one described by Borah et al. 
(1981), which simulates the supply of water and sediment runoff from 
adjacent upland areas. At present the model is restricted to 
consideration of noncohesive bed materials only. It is also restricted 
to down channel streamflow, and cannot consider the effect of transverse 
currents. 

Parts 2 and 3 of this report provide a detailed description of the 
model. Part 4 discusses the validation of the model on sets of 
laboratory and field data. Coding details are given in Addendum 3. 
Some of this material has been presented in an earlier report (Borah, 
1979). 
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2 MODEL FORMULATION 
2.1 EQUATIONS OF MOTION 

The present model treats the time-dependent problem of one- 
dimensional sediment routing in alluvial channels. The hydraulic 
functions driving the sediment movement are the local flow discharge, 
stage, and cross-sectional flow area. Therefore, an unsteady flow 
algorithm is required to compute the time and space distribution of 
those flow parameters in the channel, given information concerning chan¬ 
nel geometry, history of inflowing water discharge, and/or downstream 
stage. 

One-dimensional hydraulic routing algorithms are usually based on 
the shallow-water equations of momentum and conservation of mass for 
sediment-laden water. In natural streams load concentrations of up to 
50,000 ppm will not change the density of the mixture by more than .3%. 
Thus, density variations may be ignored. Furthermore, within that range 
of concentrations, the surface waves propagate with a velocity that is 
practically unaffected by the erodibility of the bed (Gradowczyk, 1968). 
The governing equations for water movement can thus be written 


8Q 
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( 1 ) 


aq . aA ^ ^ (2) 
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where A is the flow cross-sectional area, g is the acceleration of 

gravity, Q is the flow discharge, q is the lateral water inflow per 

w 

unit length of channel, is the bed slope, is the friction slope, t 
i.s t ..me, V is the mean flow velocity, is the velocity component of 
lateral inflow in the main flow direction, x is the horizontal distance 
along the channel, y is the flow depth, P is the momentum coefficient, 
.and p i.s the water density (Fig. 1). 

Two approximations to Eqs. 1 and 2 have found wide application in 
unsteady flow routing. They are the diffusion and kinematic wave ap¬ 
proximations (Ponce, Li, and Simons, 1978). The first assumes that the 
inertia terms are negligible, while the second assumes that both the 
inertia and pressure terms can be neglected. Both approximations have 
I'f-en shown to be satisfactory in a variety of cases. Solutions of the 
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complete sha 1 1 ow-wa t c r iNjiuit i ons and Iht'ir i iri.tl i ons have 

been extensively discussed in the literature (•'liiier anti Yevjevich, 
1975) and tliey will i,ol be considered in this pijar Anv itceplable 
hydraulic algoritlun will snflice since the waU'r inov. meru is assumed to 
be uncoupled from the si'dimenl processes. 

The third equation t t motion is given by the t (iiiserval ion of mass 
equation for sediment. As shown in Aildendiiiii i , this eijuation can be 
expressed in the form 


SQ 

cAc ^ _s 
8t 9x 


dt ^ dr 
91 91 


where c is the avenge voinriic concent rat ion, is the sediment voiume 
flux, A is the bed po ro;.:, ‘ y, v is the active wititn of the bed (i.e., 
that portion of tli-,' lu-d '',’idth in which erosion or deposition takes 
place), z is the local ‘■•ed elevation, r is the nel sediment volume flux 
through the suspended-bed load interface, and is the lateral inflow 
of sediment per unit channel length. 

In Eq. 3 the third term represents the volume raU’ of sediment 
scour (or deposition) per unit length of channel bed. The fourth term 
;is envisioned as the ti.me rate of net local exchange between the sus¬ 
pended and bed load zenes. An expression defining this term is needed 
for solving Eq, 3. For simplicity, the following [v.camelnc exchange 

equation is adapted (Whitham. 1974) 


9r 

9t 


1; Ah'! 
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c-kj(R-r )c] 


(4) 


Here and k.^ are roust.ints, T is the concentration at transport capa¬ 
city, corresponding to v/iiich, K is the net sediment exchange between the 
suspended and bed zones. Eij. 4 was so constructed in order to (i) 
incorporate in the solution some of the nonlinearity nndotibtly present 
in the exchange process, and (ii) p-reserve the hyperbolic character of 
the sediment continuity equation. Although Eqs. 3 and 4 c.an be solved 
by succes.sive iteration'■ to obtain c and r, -a simpler solution results 
from the following observation. Very near equilibrium the right hand 
side of Ecp. 4 vanishes for all practical purposes. That is 
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k^A{(T-c)r-kj(R-r)cl^ 0 . 


(5) 


At the same time T deviates slightly from its equilibrium value. Thus, 
from Eq. 5 results 

^ ^ . ( 6 ) 

3t [T+(kj-l)cl2 at 

Also this expression is assumed to hold in slowly varying flow condi¬ 
tions. Approximating by AVC over a small time interval, and combin¬ 
ing Eqs. 3 and 6 gives the sediment continuity equation used in this 
model 



9c ^ ac _ %t 

8t s 3x A(l+f ) 

s 

(7a) 

where 

V = ^ 

s (1+fg) ‘ 

(7b) 


k^RT 

^s ■ A(T+(kj-l)cl^ ’ 

(7c) 


%t = % • if 

(7d) 


Eq. 7a is a quasilinear hyperbolic equation governing the propagation of 
the sediment concentration waves. Its right hand side represents the 
lateral sediment inflow contributed by runoff, tributary channels, and 
bed scouring. The limited experience gained so far with the model 

-3 

indicates that k^ is of order 10 Thus Eq. 7c has been approximated 

by, 


f = 
s 


CEL _ , 

AT[1-0.999c/T]2 


(7e) 


where CEL is a user supplied parameter. This parameter should be ad¬ 
justed to match the time of arrival of the observed sediment load peak 
at the channel outlet. Eqs. 7b and 7e show that the celerity of the 
sediment wave, V^, is controlled by the local hydraulic conditions and 
sediment properties since capacity depends on these parameters as well. 
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Eq. 7a can be solved by means of the niethoii of cha racter isI i cs. 
From this equation and the total differential of the sediment 
concentration, the following characteristii equations result 


^ V 

dt s ’ 


(8a) 

dc _ ^St 

^st 

V ■ 

Q s 

(8b) 

dt “ A(l+f ) 
s 


equations show that 

Eq. 7a possesses only one 

system of forward 


characteristics. Accordingly, this equation cannot be used in situa¬ 
tions where there are flow reversals. Integrating Eqs. 8a and 8b with 

the initial condition c = c(x ,t ) gives 

o o o 

X 

X 

o 

"" -1 

t = t + / V (l+l )dt . (9b) 

o s 

X 

0 

These integrals are used to track the evolution of the sediment waves 
across the characteristic plane. The concentrations existing on all the 
characteristics at the time of their arrival to the downstream boundary 
define the outflow sedimentgraph. 

It is interesting to note that the preceeding equations reflect the 
expected behavior. For instance, it is a recognized fact that the 
streamwise velocity of sediment particles always lags behind the velo¬ 
city of the surrounding fluid (Francis, 1973). This velocity lag ranges 
from very small values for silt particles in suspension, to quite large 
differences in the case of coarse sands ami gravels. Thus, the celerity 
of a sediment-1 oad wave will lie smaller, in general, than the celerity 
ol the carrying flow wave. This tremi is also predicted by Eqs. 7b and 
7e, for they indicate that waves of coarse material will travel slower 
than waves of finer material given that the carrying capacity of a 
stream increases as the sediment size decreases. 

This celerity lag is strictly a function of local flow and sediment 
Iiroperties, and it should not be confused with the differences sometimes 
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observed between the time of arrival of the flow and sediment load 
peaks. In fact, the sediment load may peak ahead, in phase, or behind 
the flow peak depending on antecedent conditions, intensity of the 
event, sediment source location, season of the year, etc. (Guy, 1970). 
To illustrate this point consider a channel reach having a deposit of 
very fine sediment on the upstream portion of the reach, and receiving 
an inflow Q^j as shown in Fig. 2a. The sediment will be entrained as 
soon as reaches a sufficient intensity, and will move along charac¬ 
teristics parallel to the flow characteristics (i.e., = V, Eq. 8a). 
As the characteristics reach the downstream boundary the water and 
sediment outflows begin to rise, generally at different rates. The 
concentration increase rate is controlled by the ratio of sediment 
entrainment to flow discharge (Eq. 8b). If the supply of loose sediment 
is finite, the sediment outflow will peak at some time t^, while the 
water outflow will continue to increase and peak at a later time t.^. 
The lag between these two peaks will obviously depend on the magnitude 
of the sediment supply, inflow rates, and distance of sediment source to 
channel outlet (x = X 2 , Fig. 2a). For simplicity, this example has 
ignored backwater effects to restrict the flow movement to forward 
characteristics. Next, consider the same channel reach but with the 
source located farther upstream and a constant base flow (Fig. 2b). 
The flows and sediment loads reaching the channel inlet (x = x^) are 
generated by characteristics emanating from the upstream reach, x^ I x ^ 
Xj. Lateral runoff, q^, is assumed to begin entering the channel at t = 
t^. At this point the outlet hydrograph will begin to rise and even¬ 
tually peak at t = t^. On the other hand, the sediment characteristics, 
which enter the channel reach after t^, will reach the outlet later than 
t,^ and will finally peak at a later time t^ ^ t^. Thus, in this case 
the sediment peak lags significantly behind the flow peak as result of 
changes in the channel water inflow and sediment source location. 

Consider now a channel with sediment moving mostly as bedload. In 
this case the term 9r/9t in Eq. 3 can be neglected. Furthermore, 
Mahmood (1973) has shown that ignoring the time derivative of the spa¬ 
tial concentration has little influence on the simulation of bedload 
propagation. If in addition the period of the flow velocity field is 
very large, the surface gravity waves can be filtered out and only the 
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2. Routing of water 


sediment 










bedload propagation waves need be retained (Gradowczyk, 1968). In such 
an instance, Eqs. 1, 2, and 3 reduce to the equations governing the 
"known discharge" approximation frequently employed in modeling propaga¬ 
tion of bed transients (Gradowczyk, 1968; de Vries, 1971; Cunge and 
Perdreau, 1973; Thomas and Prashun, 1977; Ponce et al., 1979). 

The foregoing observations demonstrate that the equations presented 
in this Section incorporate the ingredients needed to simulate the 
actual physical process. These equations are complemented with a number 
of ancillary algorithms discussed in the following sections. 

In natural channels there is usually a wide gradation of sediment 
sizes. Although the finer particles comprise the bulk of the load, the 
channel evolution is controlled by the coarser particles (i.e., armor¬ 
ing). Moreover, different sizes are transported at different rates as 
pointed out earlier. It is thus important to predict the movement of 
the individual particle sizes encountered in the channel. The present 
model divides the size range in a suitable number of size fractions, and 
then uses Eqs. 7b, 7d, 7e, and 9 to track the movement of each fraction. 
2.2 ANCILLARY ALGORITHMS 

2.2.1 Sediment Transport Formulas 

These formulas are used to determine the potential carrying 
capacity oi a specific flow. Different capacities can be expected for 
different particle sizes, and not all L.isting formulas perform equally 
well for all sizes. In the present model several formulas are used that 
are framed for easy use in digital applications, and require only 
information on the hydraulic parameters of the carrying flow. They are 
the total load formula of Yang (1973) used to estimate the transport of 
very fine to coarse sands (0.1-2 mm), a duBoys-type bedload formula 
(Graf, 1971) used with fine gravel particles (2-4 mm), and the Meyer- 
Peter and Muller bedload formula (Mcyer-Peter and Muller, 1948) used 
with particles in the medium to very coarse gravel range (4-65 mm). 
These formulas are presented in Addendum 2 in the forms used in this 
model. The Yang formula was found to give very reliable estimates for 
flows carrying sands in the indicated size range, including flows 
transporting sediment mostly as bed load (Alonso et al., 1980). The 
duBoys-type formula and the Meyer-Peter and Muller formula were selected 
because they were developed from data in the specified size ranges. 
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Nevertheless, the user may replace these formujas by others he may deem 
more appropriate for a particular simuJaLion. In particular, simple 
relationships between the sediment transport rate and the flow condition 
developed from in situ field surveys should be given preference, 

2.2.2 Residual Capacity. Composition Material in Transpo rt 

Whenever the bed material consists of a mixture of different 
sediment sizes, the particle size composition of the sediment discharge 
usually differs from that of the bed. Therefore, the transport must be 
characterized by a load rate and a size distrib'ition analysis. Although 
a flow may have the potential to transport sediment, its residual 
capacity or ability to carry any additional load depends on the sediment 
material already present in the flow. Consider, for instance, a flow 
carrying a load c^ of uniform size d^ and let be the corresponding 
potential capacity of the flow. Then 


= h • - T, (C,/T,) 

is the residual capacity of the flow for that size material. The last 
term in this expression represents that portion of i^ already consumed 
by the m.aterial in transport. Similarly, if the later were of a dif¬ 
ferent size the residual capacity for the size d^ materia, , in the 

presence of the load C 2 , would oe 

T,., = T, - T, ICj/T^) , 

where envisioned as that part of T^ depleted by the load 

c^. If both sizes were simultaneously present in the flow, then 

"^rl = "^l ■ ■ "ri • 

This expression can be generalized to any size fraction, d^ , and for an 
arbitrary number, n, of load fractions c^, c^, ..., c^, as follows 

T^i = T. - t.(c^/t^) - 

- T,(c,/T.) - ••• - T.(c^/Tj, (10) 


T . 
r 1 


= T. 
1 


n 

.j = l 


(c./T.)| -AT., 
.1 J > 


j = 1 , 
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( 11 ) 
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n 

The quantity within brackets, 1-1 (c./T.), represents the portion of 

j=i ^ 

the potential capacity taken up by all the size fractions in 

transport. Therefore, the quantity within brackets, A, is the remaining 
capacity for transporting additional material of size d^. For a given 
sediment load, A is the same for all fractions. depends uniquely on 

the local flow and the properties of the d^ fraction, while depends 
on all these parameters and on the size composition of the sediment 
load. 

When A > 0, any size fraction available for entrainment at the bed 
surface, and for which 0, will be removed by the flow and added to 

the same sediment size class already in transport. Thus, A > 0 identi¬ 
fies an eroding bed condition. Similarly, when A < 0 the stream carries 
a load in excess of its potential capacity and will deposit the excess 
sediment material on the bed. Therefore, A < 0 characterizes an ag¬ 
grading bed condition. When A = 0 there is no load change and the 
transport process remains in a pseudo-equilibrium condition. By contin¬ 
uously tracking the value of A, the dependence of the individual resi¬ 
dual capacities on the composition of the sediment load and its exchange 
with the bed is readily simulated. 

The size composition of the total load is continuously up-dated by 
adjusting the concentration of individual fractions according to the 
composition of the material removed from or added to the bed. This 
procedure is explained later in the report. 

2.2.3 Bed Compo sition 

In cohesionless beds the material available for entrainment is 


essentially that exposed at the bed surface. As dunes and ripples move 
slowly downstream they continuously mix all the sediment they contain. 
The space occupied by those bed features may thus be regarded as a 
mixing zone below which the bed material remains undisturbed (Fig. 3a). 
Where the sediment contains a mixture of different sizes, the slowly 
moving coarse material tends to collect at the base of the mixing zone 
thus forming a lense of large grains. Linder some flow conditions the 
coarse material in the bed may become immobile. In this case, the flow 
washes the finer particles out of the mixing zone leaving an armor coat 
that protects the underlying material. During the scouring of the finer 
materials the exchange between bed and flow takes place in a thin layer 
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Fig. 3. Schematic representation of bed processes 









at the bed surface, identified here as the active layer. Several of 
these layers will be scoured away while the mixing zone is degrading. 
When the bed is armored, the last active layer becomes the armoring 
coat. Conversely, during the process of bed aggradation several active 
layers will be deposited on the bed forming a new mixing zone. 

To model the above processes the mixing zone is pictured as a band 
of constant thickness divided into several layers (Fig. 3b). The layer 
in contact with the flow is always referred to as the active layer. The 
thickness, porosity, and size distribution of this layer can vary 
throughout the simulation, but the layer is assumed to be homogeneous 
within itself at any given time. 

When all the material within the active layer moves, a reasonable 
upper bound of its thickness is obtained, from volumetric considera¬ 
tions, as 


ALT 


100 

P 

n 



( 12 ) 


Here d^^, P^, and are the size, percentage, and porosity of the 
coarsest fraction in the siniiment mixture. For instance, in a uniform 
bed material with A^^ = O.SO, Eq. 12 gives ALT = 2d^ which is in agree¬ 
ment with the bed load thickness proposed by Einstein (19.'>0). However, 
when the active layer becomes an armor coat, its actual thickness may be 
considerably less than that predicted by Eq. 12 because the bed may be 
armored by particles smaller than <1^^. In that rase, it is proposed to 
compute the layer Lhi'kness from 


ALT 


100 

n 

2 P 
i=2 



1-A 
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(13) 


where d^ is the smallest grain fraction the flow cannot transport (i.e., 
T.=0, i = ,n)■ Eq. 13 is also a better measure of the active 
layer thickness when some of the fractions in the bed layer cannot be 
eroded by the flow. At low discharges only the smaller fractions will 
be set in motion and Eq, 13 will thus predict a thin active layer. This 
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is intuitively correct since little bed material will be scoured by a 
low flow during a simulation time step. As the discharge increases, the 
coarser fractions, usually present in small percentages, will be en¬ 
trained and Eq. 13 will predict a thicker layer. This behavior is in 
agreement with the fact that a greater depth of bed can be sorted 
through by a higher flow in the same amount of time. Eq. 13 is thus 
adopted as the general expression for the active layer thickness. This 
equation incorporates Eq. 12, in the limit, by letting £ = n when all 
fractions are in motion. The sediment contained in the active layer is 
the only material available for erosion during a simulation time step. 
When the bed is armored no erosion can occur until the flow develops the 
necessary T^ for the smallest size fraction present in the armor coat. 
When this happens the armor coat becomes again an eroding active layer. 
If deposition of a certain amount of sediment occurs during simulation, 
this material is added to the bed and a new active layer thickness is 
computed based on the new mixture composition. 

In order to account for the time evolution of the active layer 
thickness, it is necessary to continuously track the size composition of 
this layer. This is done by introducing the following accounting algo¬ 
rithm. Consider a well mixed active layer with three size fractions 
dj<d 2 <d.^ being eroded by a flow with sufficient transport capacity to 
scour all three fractions (Fig. 4). The bed scouring is imagined to 
begin with the entrainment of the d^-particles exposed at the bed sur¬ 
face. Next, the d^-particles at the bed surface are removed, followed 
by the d^-grains hidden underneath the d^-grains. Finally, the d^-par- 
ticles are washed off the bed surface followed by the d^ and d^'parti- 
cles underneath them, and by the d^-grains uncovered by the removal of 
the d^-grains. Thus, one d^-particle is removed for every d^-particle 
entriined by the flow, and one d^ and two d^-grains are associated with 
the removal of every d.^-grain. This ordering can be readily extended to 
a y number, n, of size fractions, and it is summarized in the following 
"entrainment frequency" matrix. 
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F = [F..] = 
ij 


^i-2 .^i-3 


wliere i, j = 1,2,3, ..., n. In this matrix the diagonal elements indi¬ 

cate that every size fraction at the bed surface is scoured once. The 
off-diagonal elements in each row indicate the number of times each 

fraction d., l^jSi, becomes available for entrainment once d. is removed 
J 1 

from the bed. Alternatively, the elements of any column, j, express the 
number of times the d^ fraction is depleted due to the entrainment of a 
fraction d^. Obviously, the set of entrainment events in Eq. 14 is one 
of many possible distributions since in actuality any number of d.-par¬ 
ticles can be associated with the removal of a d^-grain. In a more 
general stochastic representation the F matrix would be replaced by the 
conditional probability matrix of the entrainment process. Neverthe¬ 
less, the proposed algorithm is adopted for its deterministic simpli¬ 
city. It should be noted that the F matrix is time invariant and de¬ 
pends only on the number of fractions used in the simulation. 

The amounts of eroded materials corresponding to the frequencies 

F,, depend on the volumes of the fractions contained in the active 
ij 

layer. To convert those frequencies to volumes let be the total 

volume of the d^-size present in the layer per unit channel length, and 

V. . the portion of V, that becomes available for entrainment when the 
ij J 

d^-fraction, occupying the volume , is eroded. Thus, 

V. , = F.. V. , 
ij iJ 1 


V. = Z F , V , 
J . ri r 
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from wtiich, 
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From this expression one oblains 
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Eq. 15 can be rewritten as 
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where, 
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are the percentages of active-layer material in each size class inter 
va 1. 


Eq. 17 defines the' element:, of a sijiiare matrix identified in here 
as the "volume enltainment matrix," v. The order of this matrix is 
equal to the total number of size classes in the active layer. Its 
diagonal elements represent the volumes of the individual fractions on 
the bed surface. The off-diagon.i I row elements contain the individual 
fractional volumes exposed bv the erosion of larger sizes. Adding up 
these elements gives the vojanie of potential erosion associated with the 

removal of the largest size in the row. On the other hand, summing all 

tiie olement.s in each coliuiui yields the total volume of each size class 
present in the active layer and available for scour (Eq. 16). These 

concepts are schematically illustrated in Fig. 5a, which depicts the 
v-matrix ab.sociaLeii with a small cluster of bed particles grouped in 
three different size classes. The shaded art'a in Fig. 5b represents all 
the material that could hr enlriiiic-d along with the third fraction. The 
shaded portion of Fig. 5c ^r-presiMUs ini.tead the total volume of the 
first fraction contained in the cluster. During bed degradation, or 

aggradation, the elements .1 the >.-m.itrix are continuously adjusted as 
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explained in the next two sections. At the end ol each simulation time 
step the total volumes of the individual Iracticuis lelt in the active 
layer are introduced into Kcj. 18 to calculate the size distribution of 
the layer. 

2.2.4 Bed Ero sion. Armoring 

Whenever A > 0 the bed is in a degrading mode. This is 
conceptually modeled by depleting the elements of the y-matrix, one row 
at a time, beginning from the smallest fraction. The diagonal elements 
are depleted first, for they are ex|)osed to the flow when simulation 
starts. Next, the off-diagonal elements are scoured one by one 
beginning from the smallest size. When the flow residual capacity is 
not sufficient to remove all the v^l'imes in a particular row, equal 
amounts (for simplicity) are removed from all the elements in the row 
until the residual capacity is satisfied. The eroded volumes are added 
to the individual fractions already in transport, and the value of A is 
recalculated. This process continues until either A is no longer 

positive definite, or the entire active layer is worked through. For 
example, the volumes *^'^21’ *^'^22’ *^'^21 ’ 

v.^l/3, '' 33 '^^’ '' 31 '^^’ '^32^^’ would be eroded, in this order, 

out of the v-matrix shown in Fig. 5. However, the extent to which these 
volumes are actually entrained depends on the degree of detachability of 
the sediment particles. These concepts are used to construct the 
following algorithm expressing the total volume of sediment eroded out 
of each size class in the active layer during a simulation time step; 


I e. , , e. . = ERO . V. ., if AT , ^ I v, , , (19a) 

iJ ij 
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r< i 
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1 1.1 fj . , ij 


AT , > ‘ V. . , 
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(19b) 
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r.) ri ij ’ 

AT ’ v , , 


( 19c ) 



i=l, 2, n. In these equations KRO is a user supplied erodibility 
parameter. This parameter governs the aetual amount of bed material 
available for erosion during a simulation time step. EKO is calibrated 
by fitting the sediment yield volume to observed data. Everytime a new 


e_ is computed the concentration of the corresponding load fraction, 
Cj, is upgraded by letting 


. . 

c“ = c . + . (20) 

1 J A 


This concentration is then entered in Eq. 10 to update A. When A be¬ 
comes nonpositive for any particular size fraction the transport is 
termed "capacity limited." On the other hand, if A remains positive 
after depleting a fraction the situation is termed "supplied limited." 
In particular, if A > 0 after considering all fractions present in the 
active layer, the bed is said to be armored and the active layer becomes 
an armor coat. fig. 6 summarizes the foregoing simulation sequence. 

After the active layer has been worked through by the flowing 
water, its volumetric composition is given by V^-Ev, i=£, £+l,...,m, 
where £ and m are the smallest and largest fractions left in the layer. 
The actual thickness of this material is then 

mV. E* 

i=£ 1 

Whenever ALT" is equal, or larger, than the thickness obtained from Eq. 

12, the later is taken as the new thickness of the active layer. In 

some instances, however, ALT^’^ turns out smaller than ALT. In these 

cases, a thickness, 6, of undisturbed material is added to the active 

layer (Fig. 7a). Because the model does not track the composition of 

the undisturbed layer, this is always formed by original bed material. 

This material contains, in general, a smaller percentage of the largest 

size class than the ALT* layer (Fig. 7b) and, therefore, 6 must be 

reduced to account for the difference. For simplicity, a simple linear 

correction is assumed given by the ratio between the percentages, and 

in the eroded and undisturbed layers (Fig. 7b). Thus, the corrected 
m 

thickness of the new active layer becomes 
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Fig. 7. Aiiiiistment of active-layer thickness 
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( 22 ) 


p‘» 

ALT = ALP'- + (5 . 

c p.i 

m 

The same criteria is used when updating an armor coat tiiickriess, except 

that in this case the ratio P^/P*'* is replaced by P^/P^. 

mm £ £ 

Finally, the local bed elevation of the channel bed is updated by 
subtracting from it the thickness of tiie eroded material, giving (Fig. 
7a) 


II E. 

new old W . , (1-A.) ' 

1=1 1 

2.2.^ Sediment D eposit i on 

When the model senses deposition (A < 0), the flow drops sediment 
on the bed during the time step. At, and the settled material is added 
to the active layer. Within the present deterministic framework, 
deposition begins with the largest sediment f.-action and continues 
through the smaller fractions until either the stream is no longer 
overloaded (A = 0), or all the fractions in transport have been 
depleted. The volume deposited out of any fraction, . cannot exceed 
its (defect) residual capacity, that is. 

Maximum D. = T = T. A . (24) 

1 ri 1 

However, whether this amount will reach the bed during the time step At 
depends on this being not less than the average time for the sediment 
particles to settle to the channel bed. Data by Jobson and Sayre (1970) 
and by Lean (1971) indicate that this settling time may be computed 
using the particle fall velocities in the quiescent fluid. Therefore, 
the actual deposition of a size fraction during the interval At is 
calculated from 
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particles fall simultaneously from the water surface. The deposited 
volume is subtracted from the material in transport to yield the 

size-class concentration 


From the preceding equations the deposition loop shown in Fig. 8 is 
constructed. 

The settled sediment is added to the active-layer fractions, and 
both materials are assumed to mix thoroughly yielding the volumetric 
composition 


V., when D. = 0, 
1 1 


I V. + D., when D. ^0, 

i = 1, 2, . . . , n. These volumes are introduced in Eq. 17 to update the 
bed composition, and this is then used in Eq. 12 to compute the new 
active layer thickness. Finally, the local bed elevation is updated by 
adding to it the thickness of the settled material yielding 
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NUMERICAL SCHEME 


To run the model in a channel reach of length L, the open domain [0 
1 X $ L] X [t ^ 0] is covered with a rectangular grid with lines paral¬ 
lel to the X- and t-axes. There is a single family of characteristics 
associated with Eq. 7a. This leads to the solution of an initial, 
one-point boundary value problem. Initial conditions are cross-section 
geometry, bed elevation, average flow velocities, flow depth, active 
width, and bed size composition. Upstream boundary conditions are time 
histories of water and sediment inflow, and size composition of sediment 
load. Let Ax and At be the length and time increments, respectively, 
separating the grid lines (Fig. 9). The coordinates of the grid nodes 

are x = mAx and t = nAt, m, n = 0, 1,2, ... The value of any var- 
m n 

iable, say A, at the node (x , t ) is designated A . Given the above 
’ m’ n ® m 

data specified along the line t = t^, the mod.i is used to compute the 
sediment discharge, load composition, bed elevation, and active layer 
composition at all nodes on the next line t = For use in the 

numerical scheme, the sediment transport equations and ancillary algo¬ 
rithms are implemented as shown in the diagram of Fig. 11. 

The volumes of the sediment fractions available for erosion during 

At are obtained from the active layer thicknesses on the line t = t as 
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1 .m 


w" 
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ALT 


100 
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1-A. 


i = 1, 2, 


n. 


(29) 


These volumes are used in Eq. 17 to compute the volume entrainment 

matrix v . 

m 

The flow routing scheme supplied by the user is used to obtain the 
hydraulic-parameter values at the nodes .1 and K (Fig. 9), and their 
averages are used to compute an average transport capacity, T. The 
concentration of each load fraction at the end of the characteristic 
path is obtained from the following discrete form of Eq. 9a 


n 

c. = c, + 

1 ,B 1 ,m 


q .Ax 
^st 


(30) 


where Q is the average of Q, and Q , and q is the piecewise uniform 

J K S L 

lateral sediment inflow over Ax. The amount of sediment entrainment, or 
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deposition, and the updated load eom t-nl ral i on, < j , are talfuiated 

using either the erosion loop (Fig. <>j or the deposit i(;n loop (Fig. 8) 
depending on the sign of A. If the character . st i i of a particular 
fraction does not pass through the no<i<* K (l ig. 9j, the i ofueiitrat ion at 
this node is obtained by linear interpolation from 
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(31b) 


In these equations Ax* and At* are given by discrete forms of Eq. 9b. A 
similar correction is applied to (Kq. 19) to determine the actual 
amount of material eroded in the interval At. The new concentrations 
are used to calculate the new volume load fractions 


n+1 _ T n+1 

'i,m+l ^i,m+l 


(32) 


and these are introduced into Eq. 18 to update the load size composi¬ 
tion. The foregoing load updating sequence is summarized in Fig. 10. 
Last, the new bed size composition and elevation are computed. When a 
pseudo-equilibrium condition is encountered (A = 0) no bed updating is 
necessary. In continuous simulation the above calculations are per¬ 
formed consecutively for each chaiitiel segment to update all variables 
over the length of the channel. The process is then repeated for each 
time interval in the simulation period. 
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MODEL TESTING 


Tests were conducted to verify the ability of the model to simulate 
various instream processes. Published laboratory and field studies were 
scanned for data suitable for testing the algorithms describing proces¬ 
ses of bed scour, armoring, and unsteady transport. Three useful sets 
of data were found. The first set was collected by Ashida and Michiue 
(1971). They performed a series of laboratory experiments to study the 
effect of sediment gradation on channel armoring. The second data set 
was collected by Lane and Carlson (19S1) in the San Luis Valley canals 
in south-central Colorado. They made measurements on the t>ed and bank 
materials that formed the canals. The beds of these canals have become 
armored over the years. The data offer an excellent opportunity for 
testing the model on a natural system. Finally, the model was tested 
using data collected by the U. S. Geological Survey on a reach of the 
East Fork River, Wyoming (Mahoney et al., 1976). These data offer the 
opportunity of checking the model algorithms on an alluvial stream under 
conditions of unsteady flow. These tests are discussed below. 

As is usual in channels with material consisting of a wide range of 
grain sizes, the bed slopes of the channels used in the above three 
studies were fairly steep. In these cases the kinematic-wave 

approximation to the equations governing unsteady flow of water is 
applicable. In this approximation the momentum equation, Eq. 1, becomes 
Q = KIN . (IT) 

where 


KIN = C^ 



1.49 j^l/6 
n 



(T4) 


is a k j iieinatic-wave parameter. In these relationships P and R are the 
wetted perimeter and hydraulic radius of the channel cross section, 
respectively, C^ is the Chezy coefficient, and ri is the Manning 
roughness factor. Eqs. ,TT and 34 are also valid when the flow is 
uniform and steady. For the purpose of simulating the East Fork River 
data Eqs. 2, 33, and 34 were solved using the kinematic-wave routing 

scheme developed by Borah et al. (1980). 

The simulations required calibration of the model parameters to 
obtain best-fit of model predictions to observations. One flow routing 
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parameter, KIN, and two sediment transport parameters, CKL and ERO, were 
available for adjustment. CEL controls the travel time of sediment 
load, and ERO gove ns the total quantity of bed material available for 
entrainment at the bed active layer. 

4.1 FLUME ARMORING STUDY 

Several laboratory experiments were carried out by Ashida and 
Michiue (1971) to study the effects of bed armoring on channel 
degradation. Various sediment mixtures were used with different mean 
diameters and geometric standard deviations. These mixtures were placed 
in a recirculating flume with a bed 2.62 ft. wjde and 6S.S ft. long. In 
each experiment a steady uniform flow was passed over the bed. The 
hydraulic conditions were chosen to purposely induce armoring. No 
additional sediment was introduced into the stream, and the rate of 
sediment collection in a trap at the end of the flume was used as a 
measure of the total sediment discharge. 

Data from Ashida and Michiue's run No. 2 were chosen to check the 
performance of the model in simulating the transport of cohesionless 
sediment mixtures and the formation of armoring. The hydraulic condi¬ 
tions used in the run are given in Table 1. In running the model, the 
flume was represented by four reaches of equal length with a point 
inflow of water at the upstream end of the channel. The initial 
particle size distribution of the bed material employed in the run is 
shown in Fig. 12. In simulating sediment transport, the graded bed 
material was represented by ten discrete particle size fractions. They 
are listed in Table 2. Fig. 12 shows the size distributions of the 
eroded material leaving the flume and of the armoring layer as reported 
by Ashida and Michiue, and the best-fit curves predicted by the model. 
The agreement between the simulated and the measured armored layer size 
distributions is good. However, some discrepancy exists between the 
curves for the eroded material. The measured curve definitely shows 
larger sizes present in the sediment load. This discrepancy may be 
attributed to inaccuracies in the estimation of sediment entrainment. 
In the transport capacity formulas an average tiactive force concept is 
implied, which excludes the effect of the large instantaneous 
fluctuations associated witli turbulent tractive forces (Alonso and 
Coleman, 1981). 
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Table 1 - Hydraulic conditions in the flume and San Lujs Valley Canal 
studies. 


Data 

Source 

Flow 

Rate 

__ 

Flow 

Depth 

Average 

Ve1o cit y 
(fl)s) 

Energy 

Slope 

Manning's 
Coefficient 

Ashida and 
Michiue (1971) 
Run 2 

1.06 

0.22 

1.81 

0.00440 

0.017 

Lane and 

Carlson (1953) 

128.0 

1 . 77 

4.00 

0.00243 

0.023 


Test Section 12 


Table 2 - Size distributions used for simulating the flume and San Luis 
Valley Canal data. 



Flume 

armoring study 

San Luis Valley 

Canal 

Size Class 

Percent 

Size class 

Percent 

Interval 

per class 

Interval 

per class 


(mm) 

interva1 

(mm) 

interval 

0.2 

- 0.3 

9.5 

0.000 - 0.149 

4.0 

0.3 

- 0.4 

9.5 

0.149 - 0.297 

5.7 

0.4 

- 0.6 

11.0 

0.297 - 0.590 

10.1 

0.6 

- 0.8 

6.0 

0 590 - 1.190 

12.0 

0.8 

- 1.0 

4.0 

1.190 - >.380 

9.2 

1.0 

- 2.0 

12.0 

2.380 - 4.760 

5.4 

2.0 

- 4.0 

18.0 

4.760 - 9.525 

13.3 

4.0 

- 6.0 

18.0 

9.525 - 19.05 

17.9 

6.0 

- 8.0 

10.0 

19.05 - 38.10 

13.2 

8.0 

- 10.0 

2.0 

38.10 - 76.20 

9.2 
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fig. 12. Size d 1 s Lr 1 bii L i on ci.r-.i;. ol'i.lined for flume armoring test 









4.2 


SAN LUIS VALLEY CANAL TESTS 


The test described in this section is based on data collected by 
Lane and Carlson (1953) on a reach of one of the canals locaied in the 
San Luis Valley of southern Colorado. These canals were constructed in 
the late 1800's in the alluvial cone deposited by the Rio Grande River 
on the floor of the San Luis Valley. Near the appex of the cone the 
deposits consist of sand, gravel, and cobbles, the size of the cobbles 
decreasing with the distance from the appex. The canals were found very 
stable, and therefore presented an unusually favorable condition for 
studies of stable canals. Several reaches ranging in length from 600 to 
2200 ft. were selected for study. Hydraulic measurements were made on 
the most regular portions of the reaches. Observations and mechanical 
analysis were also made of the bed and bank materials forming the test 
sections. Samples of the material in which the canal was constructed 
were obtained by excavating into the bank of each section. Mechanical 
analysis of the material forming the bed surface layer disclosed that in 
the stable sections the finer material had been removed from the layer 
and an armor coat of coarser material was left. 

Data from the test section No. 12 was selected to simulate the 
development of the armor coat. The hydraulic conditions observed in 
this section are presented in Table 1. The bank and bed material size 
distributions are shown in Fig. 13. For simulation purposes, a 600 ft. 
long reach divided into four equal segments was assumed, with a constant 
inflow of clear water. The initial bed-material size composition, 
assumed equal to that of the banks, wa.s represented by the distribution 
listed in Table 2. The size composition of the active layer at 2.5, 
12.5, 37.5, and 87.5 hours are shown in Fig. 13. As time increases, the 
layer is quite rapidly depleted of particles 10 mm in size or smaller, 
and the distribution in this si/.e range approaches assymptotically^ the 
measured curve. However, the simulated distribution in the 10 to 75 mm 
,ze range differs from the composition. This difference may be 
explained in part by the use of a constant rate of discharge. Although 
information related to the actual history of flows in the canal is not 
■available, Lane and Carlson (1953.) mentioned that the canal had 
sustained flows considerably alicive those observed duiing their study. 
Therefore, the observed bed-layer composition was most probably shaped 
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Fig. 13. F.volution of Llie s i zo il i s I r i but i on in to[) layer of San Luis Valley 
Cana 1. 
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by flows 

having transport 

capacitles 

larger than 

that used in 

111 is 

simulation 

Another reason 

may be the 

use of average 

1 1.11 1 i ve for c 

OH as 

indicated 

in the previous 

section. 

Nevertheless, 

there is a 

clear 

tendency 

for the simulated curve 

to approach 

the measured 

s i ze 


distribution over the entire size lange, 

4.3 EAST FORK RIVER PROJECT 

A detailed program of hydraulic and sediment transport data 
collection was undertaken by the U.S. Geological Survey during 1975 on a 
reach of the East Fork River near Boulder, Wyoming. The study reach is 
approximately 3 miles in length and has a tributary area of about 180 
square miles. About half of this area lies within the Wind River 
Mountains. The other half is provided by the area drained by the Muddy 
Creek tributary. A map of the study reach showing the position of the 
principal gaging stations is given in Fig. 14. The river in the study 
reach is about 60 ft. wide and 4 ft. deep at bank full stage. Most of 
the water during high flows comes from melting snow of the mountain 
area. In the East Fork River the bed is gravel on the riffles and bars, 
but coarse sand constitutes the bulk of the bed load. The bed material 
in the Muddy Creek is sand but much finer than the sand in the East 
Fork. Stage recording gages were installed at sections B-1, B-5M, and 
B-17. Sediment inflow to the reach was measured at sections B-1 and 
B-5M by sampling suspended sediment with a standard DH-48 hand sampler 
and determining the bedload discharge with a Helley-Smith bedload 
sampler. The bedload discharge past the section B-17 was measured with 
both a bedload trap and a Helley-Smith sampler. Leopold and Emmett 
(1976) and Mahoney et al. (1976) describe in detail the stream, the 
bedload trap, and the data (ollected during 1975. The data included 
cross-sectional surveys, flow rating curves, water temperature, sediment 
transport rates, and particle-size distributions of bed load material. 
This data set was also simulated by Bennett and Nordin (1977). Their 
input data reduction procedures were followed to some extent in the 
present test. 

In this simulation, the study reach was divided into the fourteen 
subreaches shown in Fig. 14. Surveys of cross sections were used to 
obtain relationships between area, wetted perimeter, and stage. Time- 
averaged active-bed widths were determined from the surveys by comparing 
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plots of bed elevation at tbe beginuiiiK and I'ud ot the s iimi 1 at i on pci- 
iod, and observing which parts ol the bed moved vertiially. As noted by 
Bennett and Nordin (1977) this prodecure laay overestimate the aitive 
width because not always all parts of th<* active width movi simultan¬ 
eously. However, there may he no better way, short ol continuous bed 
monitoring. The profile ot six of the cross sections, ai. measured on 
June 2, 197.5, are plotted i .i Kig. 15. Also plotted .ire the estimated 

active-bed widths, and the stages observed during June 7. As these 
figures reveal, overbank flow occurred at several sections along the 
reach during the high flow periods, including the inlet sections K-1 and 
B-5M. 

Flo’.' input to the system was provided by point loads at sei’tions 
B-1 and B-5M. These point loads were generated by converting the stage 
readings to flow rates using the rating curves constructed from the 
stages and flow discharges measured at those sections The rating 
curves reported by Mahoney et al. (1976) are shown in Figs. Ibu and I6b. 
The sum of the inflow hydrographs measured during the simulation period 
is about equal to the outflow hydrograph observed at section B-17, 
indicating that water was canied by essentially translatory waves. 
This fact lends further support to the use of a kinematic-wave routing 
scheme. 

Inflow of sediment at section.s B-1 and B-5M were obtained from 
sediment-flow rating curves. fhese wt>re constructed by relating the 
measured sediment loads to water discharges obtained from the above flow 
rating curves. Suspended sediment did not contribute signii''"antly to 
the total input load and, therefore, the sediment discharges were based 
solely on the bedload data. When plotted in log-log scales the data 
exhibited considerable scatter but displ.iyed a linear trend. Linear 
interpolation of the logarithmic values yielded the rating turves shown 
in Figs. l6c and l6d. The data scatter resulted in the low correlation 
< ■fficients indicated in the figures. A total of ten different frac¬ 
tions with mean diameters 0.12, 0.,;o. C.h. 1 .1), 2.0, u.O, 8.0, 16.0, 

)2.0, and 6^.0 mm were selected to repr< .eiit the seilimeiil i ii the bed and 
in transport. The percent ages at in.iteri.il associ.ited with each size 
were determined from the s,imple<l • i/e d ■ st i i hut i uis Jhese percentages 
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were used to convert the total bedloads to sediment transport rates for 
each particle-size fraction. Similar analysis was performed for the 
data collected at section B-17. These data are used in comparing 
measured and predicted water and sediment discharge 1ioin the study area. 

A thirty minute time step was used in modeling the East Fork River 
for a period of twenty-two days from May 29 through June 19, 197‘i. The 
first step involved calibration of the flow jiarameter to best fit the 
observed outflow hydrograph. A comparison of simulated and recorded 
flows at section B-17 is presented in Fig. 17. Agreement is satisfac¬ 
tory at low and intermediate fb-ws, but the measured values exceed the 
simulated values at high flows. The reason iur this discrepancy is the 
occurrence of overhank flows which were not accounted for in the con¬ 
struction of the flow rating curves (Figs. 16a and I6b). 

The next step after obtaining the hydraulic results was simulation 
of the sediment transport. The simulation was performed a number of 
times, adjusting each time the sediment parameters to obtain the best 
possible agreement between simulate<l and observed values. Fig. 18 shows 
predicted and recorded bedload discharges at section B-17 for the 
twenty-two day simulation period. The dashed lines joining the data 
points are imaginary lines used to approximate the shape of the measured 
sedimentgraph. The observed total bedload outflow between June 1 and 
June 19 is 1889 tons. The predicted value is 1942 tons, less than three 
percent higher. In spite of this .agreement, there is a tendency for the 
simulated bedload rates to be lower than the recorded values at high 
flows, and higher during the recessions. The differences are most 
noticeable during the first high flow period. A possible explanation 
for these discrepancies is provided by field observation in the East 
Fork River recently reported by Meade et al. (1981). Over the years 
they have observed that, during the runoff season, some sections of the 
reach are scoured while other sections are filled, resulting in a 
nonuniform storage of movable bed material along the stream. 
Consequently, the relation of hedload rate to water discharge varies 
significantly from one part of the reach to another. Meade et al. 
(1981) noted that iminediatelv liownsLream of a storage area the bedload 
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transport may increase steeply with the initial increase in water 
discharge and later, as the stored movable material becomes depleted the 
transport rate rapidly decreases, resulting in a multi-valued rating 
curve. They also mention that this behavior has been indeed observed at 
the gaging station B-17 (Leopold and Emmett, 1976). This behavior 
cannot be exactly reproduced by the present model, which assumes a 
supply of transportable material unitormly distributed along the 
streambed, and single-valued relations between bedload transport and 
water discharge at all sections. Nevertheless, the predicted 
sedimentgraph displays the rising and falling trends observed in the 
recorded bedload during the twenty-two day time span. 

Fig. 19 shows a comparison of simulated against sampled size dis¬ 
tribution of the bedload at the beginning, middle, and end of the first 
high flow period. The model predicts fairly well the distribution of 
the coarser materials, and the shift of the d^^ during the event. The 
recorded values, however, indicate less mobility of material in the 
finer size range than predicted. There could be several reasons for 
this disagreement. For one thing the kinematic-wave approximation used 
in routing the flow may be distorting the local energy slopes, resulting 
in overestimation of entrainment threshold conditions. Or, there might 
have been more fine material being carried in suspension which could not 
be sampled as part of the bedload, but which was actually accounted for 
by Yang's total load formula used in the simulation. A definite answer 
requires further investigation. 
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S CONCLUSIONS AND KF.COMMKNDAT 1ONS 

5.1 CONCLUSIONS 

1. A one-dimensional numerical model has heen developed for simulating 

the movement of well graded sediments through a stream network. 

2. Hydraulic routing is performed by using any acceptable algorithm 
supplied by th(' user because the water movement is assumed 
uncoupled from the sediment process. The model can be used in 
conjunction with any suitable sediment yield model to supply the 
water and sediment runoff from lateral areas. 

3. The sediment routing scheme' is b.ised on the physical proces.ses 

governing the mechanics of sediment movement in alluvial channels. 
The model recognizes the effect of bed and suspended load 

interaction on the total load movement, can simul.ite bed armoring, 
changes in bed elevation, and longitudinal sorting of eroded 
material. 

4. The applicability of the model is restricted to noncohesive 

materials, relatively stable channel geometries, streams with 
negligible in and out-of-bank transport, and flows in which 
transverse currents may be ignored. 

5. The model gave satisfactory results when tcste<l on laboratory data 
from a flume armoring study, and field data from the San Luis 
Valley Canal, Colorado, and the East Fork River, Wyoming. These 
test , tend to indicate that the model adequately simulates the 
transport of graded cohesionless sediments, including the effect of 
armoring. 

5.2 RECOMMENDATIONS 

i . It IS K’commended that the channel model be further tested against 
a vaiiety of real situations witli i.perial emphasis on the scour, 
deposition, and transport of noncohesive materials. 

2. It is recommended that the model be further developed and refined 
to include the following capabilities: (a) improve the one¬ 
dimensional representation by separating flows in the incised 
channel from flows over flood plains; (b) account foi' in and out- 
of-bank transport, and lateral distribution of bed-material 
properties and hydraulic conditions; (c) predict the variation of 
lateral bed slope and lateral sediment sorting around channel 
bends; and (d) simulate sediment retention by grasses and other 
vegetative covers. 


J.52 







3. The channel model should be tested using hy[iothelila1 situations to 
confirm that the model does respond in a realistic manner. for 
instance, these tests may include the folkn.ing channel-stabi1ity 
related applications: fal Consolidate the channel model with 
continuous sediment yi<*ld aiul bank-stability models. Kun the 
consolidated models for a period of a few years, and predict the 
size and grade of channel needed to maintain a bank height-slope 
that is stable for a given stratigraphic condition, (b) Run the 
consolidated model for a combination of unstable bank and steep 
grade and observe what combination of bed armoring and/or grade- 
control structures are predicted to slabili/,e the channel. (c) 
Select a range of storm events and use the consolidated model to 
study slough of h.mk material, and find channel width and/or 
armoring coat that is needed to prevent erosion of slough material. 

4. It is recommended that data gathering efforts he continued to 
provide an adequate base for further model development and 
validation. 
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ADDENDUM 1: SEDIMENT CONTINUITY EQUATION 


Consider an unsteady streaniflow carrying sediment down a channel 
with arbitrary cross-sectional geometry and alignment (Fig. 1.1). The 
volume concentration of the sediment-laden flow, c, may vary from point 
to point, and as a result of convection, from instant to instant at any 
point. The lateral volume rate of sediment inflow, q^, can vary with 
both space and time. The continuity equation for an infinitesimal unit 
volume in the neighborhood of a fixed point is 


3c 

3t 


+ 



- 0 . 


( 1 . 1 ) 


in which Uj^ is the instantaneous local velocity component of the 
sediment-laden fluid along the Xj^-direction. Repetition of the sub¬ 
script k in a term implies summation over the three orthogonal coordi¬ 
nate directions. In a turbulent flow, time averaging Eq. 1.1 yields 


3c 8 

3t 3x. 
k 


(c u^^ + C-U^ ) = 0, 


( 1 . 2 ) 


where the overbars indicate temporal means, and the primed terms repre¬ 
sent turbulent fluctuations. Let assume the mean cross sectional area, 
A, of the channel divided in two parts. One part, A^, is occupied by 
the sediment carried in suspension, and another, A^, is occupied by the 
sediment transported as bedload. It should be noted herein that the 
orientation of the coordinate system is arbitrary and, in general, A is 
a function of as well as of time. For convenience, the coordinate 
system will be chosen so that the direction normal to A coincides with 
the streamwise direction x (Fig. 1.1). The continuity equation for the 
suspended-load section is obtained by integrating Eq. 1.2 over A^, to 
obtain 


;/ 


3c 


- dAj + // (Cjii + cju’) dA 


‘I 


(1.3) 
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Using Leibnitz's rule: 


1 — _ _ 
^ SS dAj - [c g^J_ + 3 ^ JJ(CjU + c;u'J ciAj 


A1 


o 


A1 


BA, 




(1.4) 


in which is the mean perimeter bounding the mean area A^ . 
dary condition for this cross section is 


The boun- 


dA 


8A 


9A, 




(1.5) 




where q^, the volume rate of lateral sediment inflow per unit length of 
, is taken positive for inflow. Expanding Eq. 1.5 and time averaging 
gives 


aAj aAj _ aA^ 

^^1 ^ ^ "i ^ ^ ^ ax" 




a, 


Substituting this equation into Eq. 1.4 yields 


1^ SS c^dAj+ 1^ SS (CjU + cju’)dA^ = J (c^q^ + ^'q^do^ (1.7) 


1 


1 


1 


Analogous to the Reynolds closure scheme for turbulent momentum trans¬ 
fer, an eddy mass diffusivity tensor, e, is introduced such that 


c:u = -c 


3 x 


( 1 . 8 ) 


Now let 


Cl - Cj + C-J-, 


(1.9) 
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and 


u = u + u*, 


( 1 . 10 ) 


where the tildes denote spatial averages over Aj, and the asterisks 
denote deviation from the spatial averages. Substitution of Eqs. 1.8, 
1.9, and 1.10 into Eq. 1.7 gives 


3t 3x 


3_^ 

3x 


JX(c 



3x 


C"'u*)dAj + 


where 

Qj = JJc^u dA^ 
A1 


ol 


( 1 . 11 ) 


represents the volume rate of suspended sediment discharge. The first 
integral on the right hand side of Eq. 1.11 denotes the rate of longitu¬ 
dinal dispersion of sediment. The second integral can be replaced by 
the sum of the sediment inflow from runoff and tributaries, q , and the 
rate of sediment volume transfer across H-H’ (Fig. 1.1) from the bedload 
zone, denoted herein as - Sr^/Ot. Longitudinal dispersion in the sus¬ 
pended zone is negligible with respect to vertical dispersion and can be 
omitted. Thus, Eq. 1.11 becomes 


3AjC^ 


3r 

3t 

3x s 

^t 


( 1 . 12 ) 


Similarly, integrating Eq. 1.2 over A^ yields the continuity equation 
for the bedload. 


3t ^ 3x 


.f (0292 + e'q') da^, 


(1.13) 
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where is the volume rate of bedload discharge, and is the mean 
perimeter bounding Two different sources contribute to the integral 
in the right hand side of Eq. 1.1.3. One is the rate of sediment volume 
transfer across H-H' from the suspended zone, and denoted Or^/Ot. The 
other contribution comes from the rate of sediment exchange with the bed 
surface, which can be approximated by 

= -(l-X)W || , (1.14) 


where k is the bed porosity, W defines the active bed width (Fig. 1.1), 
and z is the local bed elevation above a reference datum. The negative 
sign in Eq. 1.14 accounts for the fact that the sediment flux across 
is negative when sediment settles out of the bedload zone during bed 
aggradation (3z/3t>0). Replacing the integral in Eq. 1.13 by the sum of 
3r2/9t and q^^ gives 


3A,c„ 3Q 3r . 

.. ?■ ^ + __2 = _2 . ci-A)W — 
3t 3x 3t ^ 3t 


( 1 . 1 '^) 


Adding Eqs. 1.12 and 1.1.*) results in 


3Ac ^^s . 3z 


ar. 

_ ^ 

3t 




(1.16) 


where Ac = A^c^ + A^c^, and ~ ^ ^2' term within brackets 
represents the net rate of sediment volume transfer across H-H'. This 
term is different from zero whenever a net amount of sediment passes 
from the bedload to the suspended load zone, or vice versa, as a result 
of bed scour or sediment deposition. Letting 


3t 




and dropping the tildes and overbars, Eq. 1.16 becomes 


3Ac ^*^s . .. 3z 3r _ 

(i.A) W ^ + 3^ - q^. 


(1.17) 
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in which is the total volume sediment discharge, and e is the total 
average volume concentration. Eq. 1.17 is the sediment continuity 
equation used in the present model, and equally applies to streamflows 
carrying sediment mostly in suspension or as bedload. 
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ADDENDUM 2: SUMMARY OF TRANSPORT FORMULAS 
This Addendum describes briefly the transport formulas used in the 
present model. Each formula is presented as the originator intended, 
but rewritten in terms of dimensionless parameters. Where the formula¬ 
tions required graphical solutions (i.e., determination of threshold 
conditions from Shield's curve), analytical equivalents (not shown here) 
have been worked out to facilitate their use in digital computation. 
Total load formula of Yang (1973); 

* = (v/u.^){io‘^ ‘^/S), (2.1) 

where 


<]) = 5.435 - 0.286 log (w .U57 log (u^./w) + 

(1.799 - 0.409 log (w d^^/v) - 0.314 log (u*/w)l 
log (VSq/w - V^Sq/w), 


y J 2.5/[log(u*d^Q/v) - 0.06] + 0.66, 0<Uj.d^Q/v<70, 

J^2.05, u,.d^jj/v i 70. 

DuBoys-type bedload formula (Graft, 1971): 


( 2 . 2 ) 

(2.3) 

(2.4) 


where 


♦ = a 02(1 - e/e^) 


a-x(S-l)3/V8'/V/2, 


(2.5) 

( 2 . 6 ) 


Empirical relationships for the sediment coefficient x ate given by 
Graft (1971). 

Bedload formula of Meyer-Peter and Muller (1948): 

<1) = 8 0^/^ (2.7) 

in which 

K = (V/u^)(f^^/8)^^^ (2.8) 

The friction factor associated with bed skin friction, fj^, is obtained 
from (Vanoni, 1975) 
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(f, + 2 log (k /4 R. ) = 1.14 - 


s' b 
37.40R 


2 log [1 + 

- 




. u k 

^ 1 , 3 <^^ 70 , 


u.,.k 


(fj^) ^ = 2 log (4Rj^/k^) + 1.14, —^ > 70, 


(2.9) 

( 2 . 10 ) 


where Nj^ is the flow Reynolds number, and Rj^ is the bed hy¬ 
draulic radius. 

In Eqs. 2.1, 2.5, and 2.7 4> is the dimensionless volume transport 
rate, 0^^ and are the mobility number and relative roughness based on 
the dj^ grain size, and the subscript c denotes threshold conditions. 
These parameters are defined as 



♦ = (Qg/W)/[(S-l)gd|Ql^, 

(2. 

11) 


0,, = uJ/(S-l)gdj^], 

(2, 

12) 

and 





11 

(2, 

,13) 

where S 

is the specific gravity of the bed material, and u,,. 

j s the 

bed 


shear velocity. 
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ADDENDUM 3: DESCRIPTION OF COMPUTER PROGRAM 

The computer program used in the simulation of the East Fork River 
data is listed below. The program consists of a main program and 
several subroutines. The main program inputs the required data to the 
system, calls the subroutines according to the computational scheme, and 
prints out the calculated results. Subroutine WROUT is a user supplied 
program that routes water through a channel reach for each time step. 
Subroutine SROUT performs sediment routing through a reach for each time 
step, and it calls in turn subroutines DUBOYS, MEYER, SETVEL, SHIELD, 
and YANG. These subroutines compute potential carrying capacities and 
sediment transport parameters. 

In order to minimi 2 e the need for core memory, the program was 
organized by taking advantage of the absence of backwater conditions in 
the present simulations. In the adopted organization the length of the 
entire channel reach is divided into a few segments, or channel blocks, 
and the time steps used for the whole simulation period are grouped into 
several consecutive time blocks. Then, the channel blocks are processed 
in sequence for each time block. Input data for each channel and time 
block are stored in disk files. A list of the important variables in 
the computer program is given in the following section. 

The codes requires 72,741 words on a Mod Comp Classic computer 
system. This is a 16-bit machine that uses two words for each single 
precision variable. The execution time for the East Fork River test is 
approximately 9 minutes. 

List of Fortran Variables 

Name Description Units 

ACCM(IF, IFA) Element of volume entrainment matrix ft^ 

at the start of the current time step. 

IFA indicates the material fraction 
exposed by the removal of fraction IF. 

AE Area of flow cross section at the ft^ 

start of the current time step. 
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ARMHT 

BEDELV(IC) 

BEDMAT(IF) 


BEDUP 


BEDWN 


CAP 


CC(IC, IF) 


CDEP(IC) 


CHEZY 

CHI 

CUIC, IF) 


Current thickness of active layer. 

Bed elevation of section IC at 
the end of the current time step. 
Vector used to store the volume of 
individual material fractions present 
in the bed layer. 

Change in bed elevation caused by 
deposition during the current time 
step. 

Change in bed elevation caused by 
erosion during the current time step. 
Volume concentration of individual 
material fractions at transport 
capacity. 

Volume concentration of material 
fraction IF passing through the 
section IC at the end of the current 
time step. 

Coefficient of power formula relating 
area of cross-section IC to water 
depth. 

Chezy's roughness coefficient. 
Coefficient for the DuBoys sediment 
transport formula. 

Volume concentration of material 
fraction IF passing through the 
section IC at the start of the 
current time step. 


ft 

ft 


ft^ 


ft 


ft 


ftVft^ 


ft^/ft^ 


L 

ft^/sec 


ft^'25/ 


ftVft^ 
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COEF(J) 


Coefficient of the entrainment 


CONC 

CPER(IC) 

CUP(1T, 

DARM 

DELT 

DELX 

DEPO 

DMM(IF) 

DPTH 

DTS 

EDEP(IC) 


IF) 


frequency matrix. 

Capacity coin, ent ration predicleil by 
transport formulas. 

Coefficient of power formula relating 
wetted perimeter of section IC to 
flow cross-sectional area. 

Volume concentration of inflow of 
material fraction IF at the start 
of the current time step. 

Size of smallest bed material fraction 
that becomes immobile in the active 
layer. 

Time taken by a sediment characteris¬ 
tic to travel the distance DELX. 
Channel length increment used in the 
computational grid. 

Volume of fraction IF deposited on 
the bed during the current time 
step. 

Representative size of sediment 
fraction IF. 

Water depth. 

Size of current lime step. 

E.xponent ol power formula relating 
area of cro;:s-secl ion IC to water 
depth. 


ppm 


ftVft^ 


mm 


sec 


ft 


ft^ 


mm 


ft 

sec 
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EPER(IC) 


Exponent of power formula relating 


ERS(IF) 


ERO 


G(IT, I, IF) 


GAMA 

GMES(IT, I, IF) 


GTOTC(IT) 


GTOTM(IT) 


I 

IF 


INLS 


wetted perimeter of cross-section IC 
to flow cross-sectional area. 

Volume of sediment fraction IF 
eroded from the bed during the 
current time step. 

Parameter controlling detachment of 
bed material. 

Calculated volume discharge of sedi¬ 
ment fraction IF passing through the 
downstream end of channel block I, 
at the end of the current time step IT. 
Specific weight of water. 

Measured inflow of sediment fraction 
IF to channel block I, at the start 
of the current time step IT. 

Calculated sediment discharge passing 
through the channel outlet at the 
end of the current time step IT. 
Measured sediment discharge passing 
through the channel outlet at the 
start of the current time step IT. 

Index identifying channel block. 

Index identifying sediment-size 
fraction. 

Number of channel subreaches in a 
channel block. 


ft3 


ft^/sec 


Ibs/ft^ 

Ibs/sec 


Ibs/sec 


Ibs/sec 
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IT 


Computation time step. 


ITCOM 

KIN 

GLAT(IT) 


MEYERP 


N 


NARM 


NFR 

NSEG 

PCBCdC, IF) 


PCF(IC, IF) 


PCFF(I, IC, IF) 


Number of time steps in simulation 
period. 

Kinematic-wave routing parameter. 

Lateral volume inflow of seciiment ft^/sec 

to channel block, during the current 
time step IT. 

Coefficient in the Meyer-Peter and 
Muller sediment transport formula. 

Number of time blocks the simula¬ 
tion period is divided into. 

Integer identifying the smallest 
size fraction that becomes immobile 
in the active layer. 

Number of representative size 
fractions used in the simulation. 

Number of channel blocks. 

Percentage of material fraction IF 
present in the active bed layer of 
cross section IC. 

Percent of material finer than size 
IF present in the active bed layer 
of cross section IC. 

Initial percent of material finer 
than size IF present in the active 
bed layer of cross section IC. in 
the channel block I. 
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PCW(IF) 


Percent of material in transport 
finer than size IF passing through 
the channel outlet. 


POR(IF) 

Porosity of bed material fraction IF. 

ftVft^ 

Q(IT) 

Computed water discharge at the end 

of the channel block I, and at the 

end of the current time step IT. 

ft^/sec 

QLAT(IT) 

Lateral water inflow to channel block, 

during the current time step IT. 

ft^/sec 

QMESdT, I) 

Measured water inflow to channel 

block I, at the start of the 

current time step IT. 

ft^/sec 

QUP(IT) 

Upstream water inflow to channel 

block, at the start of the 

current time step IT. 

ft^/sec 

RESCAP 

Residual transport capacity of an 

individual material fraction. 

ftVft 

- 

expressed as volume of dry sediment 

per unit length of channel. 


RHB 

Hydraulic radius. 

ft 

SCAP(IF) 

Potential transport capacity of 

material fraction IF, expressed as 

volume of dry sediment per unit 

length of channel. 

ftVft 

SLN 

Length of channel block. 

ft 

SLOPE(IC) 

Channel bed slope at section 1C. 

ft/ft 
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SNU 

Average value of the kinematic 

viscosity of water for the simu¬ 
lation period. 

f t^/sec 

SPGR 

Specific gravity of sediment material. 

- 

SUM 

Summation term in £q. 11. 

- 

TAG 

Average unit tractive force. 

Ibs/ft2 

TBM 

Total volume of bed materia] 

contained in the active layer per 

unit length of channel. 

ft^ 

TC 

Critical unit tractive force. 

lbs/ft2 

TCA(IF) 

Sum of all the elements in row IF 

of the volume entrainment matrix. 

ft 3 

TEMP 

Average water temperature for the 

simulation period. 

Fahrenhei 

TGIN(IF) 

Vector used to store the measured 

volumes of all fractions that enter 

the channel during the simulation 

period. 

ft3 

TGM(IF) 

Vector used to store the measured 

volumes of all fractions that leave 

the channel during the simulation 

period. 

ft^ 

TGMES 

Total sediment yield measured at the 

channel outlet. 

lbs 

TGO(IF) 

Vector used to store the computed 

yield of individual fractions. 

lbs 
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TGOUT 


Total sediment yield computed at 
the channel outlet. 


lbs 


UST 

VEL 

VS(IF) 

WATMAT(IF) 

WEP 

WIDTH(IC) 

XSI(IC) 


Bed shear velocity. 

Average velocity of flow. 

Settling velocity in quiescent water 
for material fraction IF. 

Vector used to store the volumes of 
all the material fractions being 
carried by the flow during the 
current time step. 

Wetted perimeter of channel cross 
section. 

Vector used to store the active bed 
width of all channel cross sections. 
Distance of cross-section IC to 
upstream boundary of channel block 
containing IC. 


ft/sec 

ft/sec 

ft/sec 

ft^ 


ft 

ft 

ft 
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C LIST OF COMPUTER PROGRAM USED IN THE EAST FORK RIVER TEST 

c innn%*********t*****************************$t*******t*** 
c 

DIMENSION TITLE(20),QMFS<264,3),QtlUT(300),GMES(264,10,3), 
«iPCFF(3,6,10 ) ,PCF<21,10) ,PCFR(2i , 1 0 > , GTOTM < 264 ) ,GT0TC(264) , 
iQS(21),INT<21),GS<21),ClI<3,6,10 >,BEDLVL(3,6),TGO<10),TGM(10), 
4TGIN<10) 

COMMON /ROUT/ I,IT,SLN,CHE2Y,DTS,ITCOM,INL,QBASE, 
iQUP<300),01(50,2),XI<50),KSI(S0).0(300,3) 

COMMON /UROUT/ SLOP,QLAT(300),QC<SO,2),XC(50),KSC(50) 

COMMON /SROUT/ SPGR,GAMA,SNU,NFR,INLS,WIDTH(10>,SLOPE(10 ) , 

6CPER(10),EPER(10),CDEP(10),EDEP(10),GLAT <300 , 1 0 ),XSI ( 10), 

4iXSC(10,10),BEDELV(10>,CI(10,10),CC(10,10),PCBI(10,10),PCBC(10,10), 
iCUP(300,10),DMM<10) ,COEF(10),VS(10),POR <10),G(264,3,10),PCU(10), 
«iERO,CHI ,CEL 
C 

C DATA INPUT 
C GENERAL INFORMATION 

READ(4,406) NSEG,DTS,ITCOM,NFR 
C PHYSICAL PROPERTIES 

READ(4,407) TEMP,GAMA,SPGR ,SNU 
READ(4,408) (DMM(IF>,IF = 1,NER ) 

C WATER AND SEDIMENT ROUTING PARAMETERS 
READ(4,409)ERO,CHI,CEL,CHEZY 
C INITIAL BED MATERIAL SIZE DISTRIBUTION 

READ<4,405) ((PCEF<1,IC,IF),IF-1,10),IC=1,3) 

READ(4,405) ((PCFF(2,IC,IF),IF=i,10),IC=1,4) 

READ(4,405) ((PCFF(3,IC,IF),IF=1,10>,IC-1,6) 

C COEFFICIENTS OF ENTRAINMENT FREQUENCY MATRIX 
COEF(l)=l.0 
C0EF(2)=1.0 
DO 280 J=3,NFR 
280 C0EF(J)^2.0*COEF(J-1 ) 

SUBU=(SPGR-1 . 0))t.-GAMA 
C POROSITY AND SETTLING VELOCITY 
DO 88 1=1,NFR 

POR (I) = l . - 0 . 245-0 , 0864/(0 . )*DMM( I ) )#)|c0.21 
D=DMM(I) 

CALL SETVEL(D,W) 

VS(I)=W/30.48 
88 CONTINUE 

BEGINNING OF TIME-BLOCK LOOP 
DO 999 N=l,4 

C MEASURED WATER DISCHARGE AT SECTIONS B-1, B-5M, AND B-17 
DO 351 1=1,3 

READ(3,333) (OMES<IT,1),IT=1,ITCOM) 

351 CONTINUE 

C MEASURED SEDIMENT DISCHARGE AT SECTIONS B-1, B-5M, AND B-17 
DO 350 11=1,3 
DO 888 IT=l,ITCOM 
QLOG=ALOGiO(QMES(IT,ll)) 

IF(II.EQ.1) EXP0NT=3.05785*QLOG-9 33462 
IFdi .EQ .2) EXPONT = l , 92638*QL0G-3.29409 
IF( II.EQ.3) EXPONT=l.06664*aL0G-5 03701 
GTOTM(IT) = 10 . OJKJttEXPONT 
888 CONTINUE 

READ(4,4ni) NINT 
URITE(S,401) NINT 

READ(4,402) (I NT (I) , (^(S (I ) , GS (I) , ( P CF (I , J ) , J = 1 ,NFR ) ,1 = 1 ,N1NT) 
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WRITE(5,402) (INT(I),QS(I),GS<I),(PCF(I,J>,J=1,NFR),I=i,NINT) 

ia=ii 

DO 309 I=1,NINT 
PCP=0.0 

DO 301 J=i,NFR 
PC=PCF<I,J>-PCP 
PCFRd , J)=PC/100.0 
PCP=PCF<I,J) 

301 CONTINUE 
309 CONTINUE 

DO 302 IF=1,NFR 
DO 303 I=1,NINT 
IFCI.GT.l) GO TO 304 
INT1=INT<1) 

DO 305 J=1,INT1 

305 GMES(J,IF,Il)=PCFR<l,IF))KGTOTM(J) 

KC=INT(1) 

GO TO 303 

304 IF<I.EQ.NINT) GO TO 306 
IB=KC-H 
IE=INT(I) 

GD=(PCFR(I,IF)-PCFR(I-l,IF))/FLOAT(IE-KC) 

DO 307 J=IB,IE 

307 GHES( J,IF,I1) = <PCFR<I-1,IF)+GD)|!FL0AT< J-KC))*GTOTM( J) 

KC=INT<I) 

GO TO 303 

306 CONTINUE 
INT2=INT(N1NT) 

DO 308 J==INT2,ITC0M 

308 GHES(J,IF,I1)=PCFR<N1NT,IF))((GT0TM< J) 

303 CONTINUE 

302 CONTINUE 
350 CONTINUE 
C 

C BEGINNING OF CHANNEL-BLOCK LOOP 
DO 201 Ial,NSEG 
INL=40 

READ<1,403) 1NLS,SLN 
WR1TE(5,403) INLS,SLN 
C GEOMETRIC INPUT FOR CHANNEL BLOCK 

READd,404) <XSI<IC),SLOPE<IC),WIDTH(IC>,CPER(IC),EPER(IC), 
iiCDEPdC) ,EDEPdC) ,IC = i ,INLS) 

WRITE <5,40 4) (XSIdC) .SLOPEdC) , WIDTHdC ), CPER (IC ), EPER dC >, 
«,CDEP ( IC), EDEP (IC) , IC= 1, INLS ) 

DO 801 IC=1,INLS 
PCP=0.0 

DO 802 IF = 1 ,NFR 
PCF(IC,IF)=PCFF(I,IC,IF) 

PCBIdC,IF)=^PCFdC,IF)-PCP 

PCP=PCF(IC,IF) 

802 CONTINUE 
801 CONTINUE 
C 

C INITIAL AND BOUNDARY CONDITIONS FOR CHANNEL BLOCK 
SLOP=0.0 

DO 11 11=^1,ITCOM 
IF(I-2)3,4,5 

3 aUP(IT)=QMESdT,l) 

GO TO 6 

4 QUP(IT)=QMES(IT,2)+Q(IT,1) 
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GO TO 6 

5 QUP<IT)=Q(IT,2) 

6 QLAT(IT)=0.0 
DO 11 IF=1,NFR 
IF(I-2)7,8,9 

7 CUP < IT , IF) =GMES (IT , IF , 1) /(QUP (IT) !t:GAMA*SPGR ) 

GO TO 10 

8 CUP(1T,IF)=(G<IT,1,IF)+GMES(IT,IF,2)/<GAMA#SPGR)>/QUP(1T) 

GO TO 10 

9 CUP(IT,IF)=G(IT,2,IF)/QUP(1T) 

10 GLAT<IT,1F)=0.0 

11 CONTINUE 

DO 13 IC-1,INLS 
IiEDELV<IC)=DEDLVl.(I ,IC) 

DO 12 IF=i,NFR 
CI<IC,IF)=CII(I,IC,IF) 

IF(CI<IC,IF) EH). 0. 0)CI(IC.IF)=CUP(1,IF> 

12 PCIiC<IC,IF)=PCHI< IC,IF) 

SLOP=SLOP+SLOPE<IC) 

13 CONTINUE 
SLOP=SLOP/FLOAT<INLS) 

Q£iASE = QUP<l) 

DO 300 IC=1 , INL 
QI(IC,l)=QUP<i) 

QI<IC,2)=QUP(1) 

XI<IC)=SLN!tcFLOAT<IC-l )/FLOAT(INI. > 

300 KSI<IC)--0 

C LIST INITIAL BED ELEUATION AND SIZE COMPOSITION 
WRITE<S,S01) I 

WRITE<5,502) <IC,<PCF(IC,IF),IF=1,NFR),BEDELU(IC),IC=i,INLS) 
C 

C ROUTE UATER AND SPZDIMENT THROUGH CHANNEL BLOCK DURING 
C CURRENT TIME STEP 

DO 101 IT=l,ITCOM 
C 

C WATER ROUTING 
CALL UR0UT2 
C 

C SEDIMENT ROUTING 
CALL SR0UT2 
C 

IF(INL.EQ.O) GO TO 129 

C RESET INITIAL CONDITIONS FOR FLOW CALCULATIONS 
DO 132 J=1,INL 
XI(J>=XC(J) 

QI(J,1)=QC(J,1) 

QK J,2)=QC< J,2) 

KSI<J)=KSC<J) 

132 CONTINUE 
129 CONTINUE 

C RESET INITIAL CONDITIONS FOR SEDIMENT CALCULATIONS 
DO 270 IF=1,NFR 
DO 271 IC = 1 ,INLS 
IFdC . EQ. 1) GO TO 272 
CI<IC,IF)=CC<IC,IF) 

GO TO 271 

272 CI(IC,IF)=CUP<IT,IF) 

271 CONTINUE 
270 CONTINUE 

IF(IT.NE.170) GO TO 101 
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WRITE(S,S03) IT 
DO 900 IC=1,INLS 
PCF<IC,i)=PCBC<IC,l) 

DO 901 IF=2,NrR 

901 PCF(IC,IF)=PCF(IC,IF-i)+PCBC(IC,IF) 

900 CONTINUE 

C LIST BED ELEUATION AND SIZE COMPOSITION AT TIME STEP 170 

WRITE(5,S0a) (IC, (PCFdC, IF) ,IF=1 ,NFR> , BEDEL V( IC> , IC = 1 ,INLS) 
101 CONTINUE 
C DISCRETIZED OUTFLOW 
QP=QBASE 

DO 315 IT=l,ITCOM 
QCC=Q<IT,1) 

Q(lT,I)=(QP+QCC)/2.0 
QP=QCC 

315 CONTINUE 
C 

URITE<5,504) 

DO 911 lC=i,INLS 
PCF<IC,l)=PCBC(IC,i) 

DO 912 1F=2,NFR 

912 PCF<IC,IF)=PCF<IC,IF-1)+PCBC(IC,IF) 

911 CONTINUE 

C LIST FINAL BED ELEVATION AND SIZE COMPOSITION 

WRITE<5,502) (1C, <PCFdC,IF),IF=1,NFR),BEDELU(IC),IC=1,INLS) 
C 

DO 913 IC=1,INLS 
BEDLVLd,IC)=BEDELVdC) 

DO 914 IF=1,NFR 
PCFFd,IC,IF)=PCFdC, IF) 

CII (1 ,IC,IF)=CldC,IF) 

914 CONTINUE 

913 CONTINUE 
201 CONTINUE 

C END OF CHANNEL-BLOCK LOOP 
REWIND 1 

C COMPUTE HYDRQGRAPH, SEDIMENTGRAPH, AND SEDIMENT YIELD AT 
C SECTION B-17 
TGOUT=0.0 
TGMES=0,0 
DO 299 IF=i,NFR 
TGOdF)=0 0 
TGM(IF)=0.0 
TGINdF) = 0 . 0 
299 CONTINUE 

WRITE<5,505) 

DO 41 IT=l,ITCOM 
QOUT<IT)=QdT,NSEG) 

GTOTCdT)=0 . 0 
DO 42 IF=1,NFR 

GdT,NSEG, IF)=GdT,NSEG,IF)*GAMA#SPGR 
GTOTCdT)=GTOTCdT)+GdT,NSEG,lF) 

TGO(IF)=TGOdF)+GdT ,NSEG,IF)*DTS 
TCMdF)=TGM<IF)+GMESdT,IF,3)JltDT5 

TGINdF)=TGINdF) + <GMESdT ,IF,i )+GMESdT ,IF,2> )!|tDTS 
42 CONTINUE 

TGOUT=TGOUT+GTOTC(IT)*DTS 
TGMES=TGMES+GTOTMdT ) )t:DTS 

C COMPUTE AND LIST CHANGES IN SIZE COMPOSITION OF SEDIMENT 
C LOAD AT B-l7 DURING CURRENT TIME BLOCK 





990 


991 


992 


43 


902 

41 

C 


999 

C 

C 

333 

401 

402 

403 

404 

405 

406 

407 

408 

409 

501 

502 

503 

504 

505 

506 

509 

510 


511 


IF(N.NE.l) GO TO 990 

IF<IT.EQ.213 OR.IT,EQ.214.0R.it.EQ.264) GO TO 43 
GO TO 41 

IF(N.NE.2) GO TO 991 

IF<IT.EQ.3 OR IT.EQ.4 OR.IT EQ.54 OR.IT EQ.99) GO TO 43 
IF(IT.EQ.100 OR IT.EQ.144.OR.IT EQ.148.OR,IT.EQ.149) GO TO 43 
IF(IT.EQ.1B9 OR IT EQ.191 OR.IT EQ.193.OR.IT.EQ.239) GO TO 43 
IF(IT.EQ.240 OR.IT.EQ.241.OR.IT.EQ.244.OR . IT.EQ.247) GO TO 43 
GO TO 41 

IF<N.NE.3) CO TO 992 

IF(IT .EQ . 21 OR . IT . EQ 22.0R . IT .EQ.24 OR . IT .E.Q.69) GO TO 43 
IF(IT.EQ.70.OR IT.EQ.72 OR.IT.EQ.119.OR.IT.EQ.120) GO TO 43 
ir<IT.EQ.165 OR IT EQ.167.OR.IT.EQ.212 OR.IT.EQ.213) GO TO 43 
IF<IT.EQ.262.OR.IT EQ.263 OR.IT EQ.264) GO TO 43 


GO TO 41 

IF<IT.EQ 57,OR.IT EQ 100 OR.IT EQ 144.OR,IT EQ 185) GO TO 43 
IF(IT.EQ.186.OR IT.EQ 192.OR.it EQ.237.OR.it,EQ.238) GO TO 43 
GO TO 41 
CONTINUE 

PCU(1)=G<IT,NSEG,1)/GTOTC(IT)*100 .0 
DO 902 IF=2,NFR 

PCW<IF)=PCU(IF-i)+G(IT,NSEG,ir)/GTOTC<IT)*100.0 
WRITE(5,506) N,IT,<PCW(IF) ,IF='1,NFR) 

CONTINUE 

LIST SEDIMENT YIELD, HYDROGRAPH, AND SEDIMENTGRAPH AT 8-17 
URITE<5,509) TC;MES,TG0UT 
WRITE<5,510) 

URITE<5,511) <IF,TGO(IF),TGM(IF),TGIN(IF),IF=1,NFR) 
URITE<S,S12) 

WRITE(S,513) (I,<aMES(I,J),J=l,3), QOUT(I),GTDTM<I),GTOTC(I), 
4.1 = 1,ITCOM) 

CONTINUE 

END OF TIME-BLOCK LOOP 


FORMAT!10F7,2) 

FORMAT<I4) 

FORMAT(14,F6,1,F8.3,10F6.2) 

FORMAT!I4,F10.2) 

F0RMAT!F7.1,F3.5,F7.3,F7.3,3F6.3) 

FORMAT!10F7,2) 

FORMAT!14,F7.1,14,14) 

FORMAT(3F7.2,F10.7) 

FORMAT!10F7.2) 

FORMAT!4F10,5) 

F0RMAT!//i5X,'INITIAL SIZE DISTRIBUTION OF BED MATERIAL AND', 

4.' BED ELEUATION FOR EACH SECTION OF CHANNEL BLOCK',14/) 

FORMAT!15X,14,10F8.3,Fi5.7 , ' FEET') 

F0RMAT!/15X,'SIZE DISTRIBUTION OF BED MATERIAL AND BED ELE'/.ATION 
4.' FOR EACH SECTION OF CHANNEL BLOCK AT TIME. STEP',15/) 
F0RMAT!/15X,'FINAL SIZE DISTRIBUTION OF BED MATERIAL AND BED', 

4' ELEVATION FOR EACH SECTION OF CHANNEL BLOCK'/) 

F0RMAT!//15X,'VARIATION OF LOAD SIZE DISTRIBUTION AT B-17 ', 
4'DURING CURRENT TIME BLOCK'/) 

F0RMAT!5X,'TIME BLOCK',12,3X,'TIME STEP',14,5X,10F8.3) 

FORMAT!//iSX,'TOTAL MEASURED SEDIMENT YIELD =',F20.5,' LBS.'/ISX 
4'TOTAL COMPUTED SEDIMENT YIELD =',F20 5,' LBS,'//) 

F0RMAT!//13X,'SEDIMENT COMPUTED MEASURED 

4,' MEASURED'/13X,'FRACTION YIELD!LBS) YIELD', 

4'(LBS) INFLOU!LBS)'/) 

FORMAT!15X,I3,F2n.S,F20.5,F20 5) 
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512 F0R«AT<//54X,'MEAS.'/16X,'TIME COMPUTED FLOW(CFS) AT 

4,'FLOW SED. LOAD<L.EtS) AT B-17'/16X, ' STEP B-i B-SM', 

4' B-i7 B-17 MEAS. COMP.'//) 

513 FORMAT(15X,I4,6F10.3) 

520 FORMAT(2X,F20,5/2X,E15.7) 

END FILE 5 
STOP 
END 
C 

SUBROUTINE WR0UT2 
C 

C THIS SUBROUTINE ROUTES WATER THROUGH A CHANNEL BLOCK USING THE 
C KINEMATIC WAVE SCHEME DEVELOPED BY BORAH £T AL<i980). A LIST OF THE 
C IMPORTANT VARIABLES USED IN THIS SUBROUTINE IS GIVEN BELOW: 

C - 

C NAME DESCRIPTION UNITS 

C - 

C IC 

C 

C INL 
C 

C IS 

C 
C 

C KSCdC) 

C 
C 
C 

C KSKIC) 

C 
C 
C 

C Q(IT,1) 

C 
C 

C QC(IC,L) 

C 
C 

C QI(IC,L) 

C 
C 

C XCdC) 

C 
C 

C XldC) 

c 
c 

c - 

c 

COMMON /ROUT/ I,IT,SLN,CHEZY,DTS,ITCOM,INL,QBASE, 

4QUP(300) ,QI(50,2) ,XI<S0) , KSI (50 ) , Q<3() 0,3 ) 

COMMON /WROUT/ SLOP,QLAT<300),QC(50,2),XC(50>,KSC(50) 

EXP---i .5 
BET=1.0/EXP 
EXPl=EXP+i,0 
EXMi=EXP-l.0 
KIN=CHEZY*SLOP**0.5 
TERM=EXP*KIN*DTS 
QL-QLATdT) 

QU=QUPdT) 


PARAMETER USED TO LABEL INDIVIDUAL WAVES. 

NUMBER OF NODE POINTS IN THE CHANNEL BLOCK 

COUNTS NUMBER OF SHOCK WAVES OCCURRING AT THE END 
OF THE CURRENT TIME STEP. 

FLAG USED TO CHARACTERIZE THE WAVE IC OCCURRING IN 
THE CHANNEL BLOCK, AT THE END OF THE CURRENT TIME 
STEP. KSC=0 FOR CHARACTERISTICS, KSOO FOR SHOCKS. 

FLAG USED TO CHARACTERIZE THE WAVE IC OCCURRING IN 
THE CHANNEL BLOCK, AT THE START OF THE CURRENT TIME 
STEP. KSI=0 FOR CHARACTERISTICS, KSDO FOR SHOCKS. 

COMPUTED WATER OUTFLOW FROM THE CHANNEL BLOCK I , AT CFS 
THE END OF THE CURRENT TIME STEP IT. 

FLOW DISCHARGE AHEAD <L=1) OR BEHIND <L=2) OF THE CFS 
SHOCK IC, AT THE END OF THE CURRENT TIME STEP. 

FLOW DISCHARGE AHEAD (L=l) OR BEHIND (L=2) OF THE CFS 
SHOCK IC, AT THE START OF THE CURRENT TIME STEP. 

DISTANCE OF WAVE FRONT 1C TO UPSTREAM END OF THE FT 

CHANNEL BLOCK, AT THE END OF THE CURRENT TIME STEP. 

'DISTANCE OF WAVE FRONT IC TO UPSTREAM END OF THE FT 

wHANNEL BLOCK, AT THE START OF THE CURRENT TIME 
STEP . 








IF(OL.EQ.0.0 AND.QU FQ 0.0) GO TO 130 
C PROJECT ALL CHARACTERISTICS TO NEW TIME LEOEL 
AC=(QU/KIN)**{<ET + QL»DTS/2 0 
QC(i,l)=KIN*AC)f.'ltEXP 
IF(QL.EQ.0.0) GO TO 102 
XC(i)=(QC<i,i)-QU)/QL 
GO TO 103 

102 XC< i )=TERM>KAC)H*EXMl/2.0 

103 ICST=1 
KSC(1>=0 
GO TO 131 

130 ICST=0 

131 IFdNL.EO.O) GO TO 150 
IS=0 

DO 104 IC=1,INL 
IA=IC+1CST-IS 

IF(KSI<IC) EQ 0 OR QI(IC,1T,GE.QI(IC,2)> GO TO 105 
C PROPAGATION OF SHOCK WAVE 
AA=(QI(1C,1)/'KIN)**I<ET 
AEi= < QI ( IC , 2 )/K IN ) ^XtBET 
AAF=AA+QL*PTS 
ABF=AB+QL*DTS 
QC< IA,1 )=KIN#AAF)t;’t:EXP 
QC< IA,2)=KIN*ATiF)K#EXP 
IF(QL.EQ.O.0) GO TO 108 
PR0D=ALP/<EXP1*<AB-AA>J(<QL> 

XC< IA)=XI (IC)+PR01)»<AEtF*i«EXPl-AAF)K«EXPl-AB**EXPl + AA«HcEXPi) 
GO TO 107 

108 XC<IA)=XI < IC) + <ai< IC,2)-Gil( IC,1 ) )#I)TS/( AB-AA) 

GO TO 107 

C PROPAGATION OF CHARACTERISTIC WAVE 

105 KSI<IC)=0 

AC=(Q1 <IC , 1 )/KjN))i<*F<£T + OL*I)TS 
QC(1A,1)=KIN#AC**EXP 
IF<QL.EQ.0.0) GO TO 106 
XC<IA>=XI{IC) + (QC(IA,D - QI(IC,1))/QL 
GO TO 107 

106 XC(IA)=XI(IC)+TERM»AC**EXM1 
C CHECK FOR NEW SHOCK FORMATION 

107 IF<KSI(IC) GT,0) GO TO 109 
IFdA.EQ.l) GO TO 140 
IF<XC(IA).LE.XC<IA-1)) GO TO 110 

C NO SHOCK IS FORMED 
140 KSC(IA)=0 
GO TO 104 

C SHOCK IS FORMED 
110 IA=IA-1 

IF<KSC( lA.) GT . 0 ) GO TO 112 
C SHOCK IS FORMED BY TWO CHARACTERISTIC WAVES 
XC<IA) = (XC(1A)+XC(lA+1))/2.0 
QC<IA,1)=QC(IA + 1 . 1 ) 

QC<IA,2)=QC(IA,1) 

IS=IS+1 
KSC(lA)=i 
GO TO 104 

C THE CHARACTERISTIC WAVE JOINS THE SHOCK AHEAD OF THE FRONT 
112 XC<IA)=XC(lA) 

QG(IA,i)^QC<IA4l,l) 

IS=IS+i 
KSC<IA)=11 
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GD TO 104 

C CHECK IF THE PROPAGATING SHOCK IS INTERSECTED HY ANY OTHER SHOCK 
C OR CHARACTERISTIC WAWE 
109 IFdA.EQ. 1 > GO TO 141 

IF(XC< lA) . L.E . XC( IA-1) ) CO TO 113 
C THE PROPAGATING SHOCK IS NOT INTERSECTED 
141 KSC<IA)=1 
GO TO 104 

C THE PROPAGATING SHOCK IS INTERSECTED 

113 IA=IA-1 

IF(KSC(lA).GT.0) GO TO 114 

C THE PROPAGATING SHOCK IS JOINED BY A CHARACTERISTIC UIAOE 
C BEHIND THE FRONT 

XC<IA)=XC<IA+1) 

QC<lA ,2)=QC(lA,1) 

QC(lA,1)= QC<IA+1,1) 

IS^IS+1 
KSC( IA>=^ 1 1 
GO TO 104 

C NEW SHOCK IS FORHED BY TWO INTERSECTING SHOCKS 

114 XC(IA)=<XC<IA)+XC<IA+1))/2 0 
QC<IA,l)=QC<IA+i,n 
QC(IA,2)=QC(IA,2) 

IS=IS+1 
KSC<IA)=2 
104 CONTINUE 

C COMPUTE OUTFLOW FROM CHANNEL BLOCK DURING CURRENT TIME STEP 
150 CONTINUE 
IAD=0 

IF< I.NL .EQ . 0 ) IA = ICST 
IFIIA.Etl.O) GO TO 122 
XEt=0 0 
QB=QU 

DO IIS J = 1 , lA 

IF<XC(IA).GE.SLN) GO TO 117 
XB---XC(IA) 

QB=QC(IA,i) 

IF<KSC<IA) .GT.O) QB=(QC(IA,l)+QC(IA,2>)/2.0 

117 IC=IA-J+1 

IF<XC(IC).LT.SLN) GO TO 115 

IAD=IAD+1 

XA^XCCIC) 

QA=QC(IC,1) 

IFIKSCdC) GT.O) QA=<QCdC,l)+QCdC,2))/2.0 

115 CONTINUE 
IFdAD.GT. 0) GO TO 118 
AI^ <QBASE/KIN>!t:*BET 

IFdT .GT , 1) AI=(QdT-l , I)/K IN) iKJKBET 

AC=AI+QL*DTS 

QA=KIN*.AC)K*EXP 

XA=TERM*AC!K#EXM1+SLN 

IF <QL . GT . 0 . 0) XA=<OA-OdT-l , I) )/QL+SLN 

118 QdT, D-QB+CQA-qBJKSLN-XBI/IXA-XB) 

INL=IA-IAD 

GO TO 123 

122 Q<IT,l)-0.0 
INL = 0 

123 RETURN 
END 

C 
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SUBROUTINE SR0UT2 


THIS SUBROUTINE PERFORMS SEDIMENT CAt.CULATIONS AT AI.L NODE POINTS 
IN THE CHANNEL BLOCK DURING CURRENT TIME STEP 

DIMENSION WATMAKi0),BEDMAT(10),SCAP(10>,ERS(10) ,ACCM<iO,10 >, 

4>TCA( 10 ) 

COMMON /ROUT/ I,IT,SEN,CHEZY,DTS,ITCOM,INL,QBASE, 
iQUP(30 0 ) ,QKSO,2),Xl( 50 >.KSI(5(l),0(300,3) 

COMMON /SROUT/ SPGR , GAMA , SNU , NE'R , 2NL5,WTDTH(1 0) ,SLOPE(10) , 
iCPER <10),EPER <10),CDEP(10),EDEP <10 >,GLAT(300,i0),XSI(10 ) , 
iXSCd 0,10) ,BEDELU(10) ,CI(10,10),CC<10,10),PCBI(10,10),PCBC<10, 10) , 
iCUP(300,10),DMM<10),COEF <10),OS(10),POR(10 >,G(26A,3,10),PCU( 10), 
4ER0,CHI,CEL 

DO 204 IC=1,INLG 

INTERPOLATE FLOW VALUES AT THE NODE POINTS 
IFdC EQ. 1) GO TO 205 
IFCXSKIC) LT XKD) GO TO 206 
DO 207 J = 1 ,100 

IF(XSI(IC) .GE.XI<J ) AND.X51(IC) .LT.XI(J+1) ) GO TO 208 
GO TO 207 

208 Qi=QI<J,l) 

IF(KSI<J) GT.O) Qi=((5I(J,l)+QI<J,2))/2.0 
Q2=QI<J^1,1) 

IF(KSI<J+1),GT.0) Q2=(QI(J+1,l)+QI(J+2,2>)/2,0 
QP = Ql + <(52-Q|i ):K<XS1 <IC)-XI ( J) )/(XI (J + 1 )-XI ( J) ) 

GO TO 209 
207 CONTINUE 
GO TO 209 

206 IF<IT.E(3. 1) GO TO 210 
Q1=QUP(lT-1) 

GO TO 211 

210 Q1=(3I<IC,1) 

211 Q2=QI<1,1) 

IF(KSI ( 1) . GT . 0 ) 02== (Q1 (1 ,1 )+QI ( 1 ,2 ) )/2.0 
QP-Ql + <(32-Ql )»XSI<IC)/XI<1) 

GO TO 209 

205 IFdT.EQ.l) GO TO 212 
G)P = QUPdT-i ) 

GO TO 209 

212 CiP = QIdC,l) 

209 CONTINUE 
C 

r UPDATE THE HYDRAULIC PARAMETERS AT THE NODE POINTS 
SLP=SLOPE(1C) 

WDTH^UIDTHdC) 

CPR-CPER (TO 
EPKEEPER(IC) 

CDP=CDEP(IC) 

E1)P-=EDEP ( JO 
EXP=1.5 
BLT = 1 0/EXP 
KIN=CHF.ZY*SL P>t;*0,5 
QE-QP 

AE-(QE/KIN)**DET 

IF(AE,LT.1.OE-S) GO TO 213 

VLL=QE/AE 

DPTH=CDP*AE!(!)i<f DP 
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UEP=CPR*AF£*)(iEPR 

RHB=AE/UEP 

HYR=DPTH 

UST=^SQRT ( 32.2*HYR*SLP ) 

IFdC.EQ. INLS) GO TO 2Ai 
DELX=XSI(IC+1)-XSI(IC) 

GO TO 242 

241 DELX=SLN-XSI(INLS) 

242 CONTINUE 
C 

C COMPUTE TRANSPORT CAPACITIES AND TERM DELTA (EQ. ii) 

SUM=0.0 

DO 220 IF=1,NFR 
GL=GLAT(IT,IF) 

CP=CI(IC,IF) 

D=DMM(ir) 

DFT=DMM(IF)/304.8 
U=VS<IF) 

IF(DMM(IF).LE.4 0) GO TO 860 
IF(DMM(IF).LT.S.O)GO TO 881 
MEYERP=a.0 

CALL MEYER(MEYERP,UST,RHU,DFT,GAMA,SNU,SPGR,WDTH,QE,UEL,CONC) 
SCAP<IF)=CONC*AE/l.0E6/SPGR 
IFCDMMdF) GT.30.0) SCAPdF) = 0.0 
GO TO 862 

860 CALL YANGCDFT,UST,SNU,VEL,SLP,U,SPGR,CONC) 

SCAP(IF)=CONC*AE/l.0E6/SPCR 

GO TO 862 

861 CALL DUtiOYS(D,DrT ,AE,UDTH,SLP ,SPGR , SNU, GAMA , CHI ,UST ,SFLOU) 
SCAPdF)=SFLOW/VEL 

862 WATMAT < IF ) ==CP*AE + GL*DELX/VEL 
IF (SCAP (IF) EG). 0 0) GO TO 220 
SUM=SUM+UATMAT(IF)/SCAP(IF) 

220 CONTINUE 
C 

C COMPUTE ACTIVE (OR ARMOR) LAYER THICKNESS 
PARM-0 0 
NARM=NFR 
DARM=DMM(NFR) 

DO 810 IF= 1 ,NFR 
IF(SCAP(IF) GT . 0 . 0) GO TO 810 
IFdF EQ . 1 ) GO TO 811 
IF(SCAP(IF-1).G1 0.0) GO TO 811 
GO TO 812 

811 NARM==IF 
DARM=DMM(IF) 

812 PARM=PARM+PCBC(IC,IF) 

810 CONTINUE 

IF(PARM.CT 0.0) GO TO 814 

DO 815 IFF=1,NFR 

IFA=NFR-lFFfl 

NARM=IFA 

DARM=DMM(IFA) 

PARM=PCPC(IC,IFA) 

IF(PARM.CT.0 . 0> GO TO 814 
815 CONIlNUt; 

814 ARMHT=10(1 ni«DARM/PARM/3n4 8 

C COMPUTE VOLUME OF BED MATERIAL FRACTIONS IN ACTIVE LAYER 
DO 813 IF ==1 ,NFR 

BEDMA!(IF)-ARMHT*PCBC(IC,IF)*UDTH/100.0 


813 





IF(l.O-SUH) 322,223,224 
CONTINUE 

DEPOSITION LOOP 
TBM=-0 . 0 
BEDUP=0.0 
DO 230 J=1,NFR 
IF=NFR-J+i 
CP=CI(IC,IF) 

IF<SCAP<IF).EQ 0 0) CO TD 231 
RESCAP--=SCAP < IF )*< 1 Q~SUH) 

DEPQ=0 0 

IF(RESCAP GE 0 0) GO TO 232 

IF<AbSCRESCAP).GT UATMAI(IF)) GO TO 231 

DEPO=ABS(RESCAP) 

GO TO 232 
DEPO^-WATHAT(IF) 

BT=^2 0!t:VS< IF )*DTS/RHb 
IF(PT.L T 1 0)DFPO-BI*DEPO 
UATMAK IF )=UATMAT< IF) - DEPCi 
AD1) = 0 . 0 

IF (SCAP ( IF ) GT 0 0 ) ADD 3>EP0/SCAf‘ ( IF ) 
SUM--SUM-ADD 

UPDATE LOAD AT NODE POINTS 
CAP-SCAP(JF)/AE 
CCC=UATMAT<IF)/AE 
XKP=CAP-0 999*02 
XKf: = CAP-0.999*GC.C 
ADD-0 . 0 

IF (CAP GT 0 . 0 ) ADD^CLL*f.AI'/ ( AF.»XKP*XKC > 

IF<AI)D LT 0.0) ADD^O 0 
FNLR=1 0+ADD 
DELT=DELX*FNLR/VEL 
IFTDELT LT.DTS) GO TO 2200 
IFdC .EQ. INLS) CO TO 233 
CCl^CI<lC+1,IF) 

CC< IC+1 , IF) = CCI+<CCC-CC1 ) *i.‘T5/DrL T 
GO TO 234 

SEDIMENT OUTFLOW FROM C.HANNFI. PLOFF; 

IFdT EQ 1) GO TO 230 
CCI=G< IT-1 , I , IF )/OdT-l , I ) 

GO TO 236 
CCI--CI < IC , IF) 

G< IT , I , IF ) -((',01 + (ECL-CCI ) *I)TG/ Df L T )*Q( IT , I ) 

GO TO 234 

0 XDEI,^1)EET*VE1./1 Nl R 

IFdC t Q . 1 ) GO TO 2201 
CCD = CCdG , IF) 

GO TO 2202 

1 CCB-CUPdT ,IF) 

2 CC( IC + 1 ,IF) + CCB+ (GCC-CCFO*DEL X/XI'El 

IFdC EQ INLS) GdT , I , IF)-C:C( JG + ) , IF )*QdT , I ) 
BFDMAT (IF ) = BEDMAT (IF > +DEPri)(cDTS/DELT 
TBM = TBM + BEDMA7 <IF) 

BEDUP-^BEDUr’ + DEPn*DTS/DCLT/UDTH/POR ( IF > 

CONTINUE 

CHANGE IN BED ELF.UATION AND Br:D LAYER COMPOSITION 
BEDEL V( IC) =Br:D[ LVdC)+BE DUP 
DO 239 IF=i ,NFP 

PCBC< IC, IF)=BEI)MAT< IF)/n'M*i00 . 0 










GO TO 204 

223 CONTINUE 
C 

C EQUILIBRIUM LOOP 

DO 240 IF=i,NFR 
CP=CI(IC,IF) 

CCC=UATMAT(IF)/AE 
CAP =SCAPaF)/AE 
XKP = CAP-0.999=t!CP 
XKC=CAP-0.999*CCC 
ADD=0.0 

IF<CAP ,GT. 0 . 0) ADD = CEL*CAP/<AE)«XKP*XKC> 
IF(ADI).LT 0 0) ADD-0 0 
FNLR = 1 0+AI)D 
DELT^ DELX*FNLR/yEL 
C UPDATE IOAD AT NODE POINTS 

IF<DELT LT DTS) CO TO 2300 
1F< IC EQ INLS) GO TU 243 
CCI = C1(IC+1 , IF) 

cc< ic+1, IF )^cci + <ccc-c:ci )*dts/dclt 

GO TO 240 

C SEDIMENT OUTFLOW FROM CHANNEL BLOCK 
243 IF(IT LQ 1) GO TO 244 

CCI=G(IT-1,I,IF)/Q<IT-i ,I ) 

GO TO 245 

2^4 CCI=CI(IC,IF) 

245 G(IT,1,ir)=(CCl+(CCC-CCI)*DTS/DELT)*Q(lT,I) 

GO TO 240 

2300 XDEL=DLLT*UEL/FNLR 
IF<IC EQ.1) GO TO 2301 
CCB=CC<IC,IF) 

GO TO 2302 

2301 CCB=-CUP < IT , IF ) 

2302 CC< IC+1 , IF)=CCB-t ( CCC-CCB ) *DELX/XDEL 
IFCIC.EQ INLS) G<IT,I,IF)=CC(IC+i,IF>#Q(IT,I) 

240 CONTINUE 
GO TO 204 

224 CONTINUE 
C 

C EROSION AND ARMORING CALCULATIONS 
C 

C COMPUTE VOLUME. ENTRAINMENT MATRIX 
DO 201 IFA=-1,NFR 
DEN=0 0 

LOOP^NFR-IF'A+1 
DO 2B2 J=1.LOOP 
IF^J+IFA-1 

282 DEN=^DFN + COEF < J )*HrDMAT ( IF ) 

DO 283 J=l,LOOP 

IF='J + IFA-1 

IF(DEN EQ 0 0) CO TO 030 

ACCM( IF , IFA)=-COEF( J)*BEDMAT( IF)/DEN*hEDMAT(IFA) 
GO TO 283 

830 ACCM(IF,IFA)=0.0 

283 CONTINUE 
281 CONTINUE 
C 

C EROSION LOOP 

DO 270 IF=1,NFR 
TCA<IF)=0.0 
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DO 279 IFA=1,IF 

279 TCA(IF)^TCA<IF>+ACCh(IF,IFA) 

278 CONTINUE 

DO 237 IF=1,MFR 
237 ERS(IF)=0.0 

DO 284 IF=1,NFR 

IF<SCAP(IF),EQ 0 0) GO TO 284 

RESCAP=SCAP ( IF )*( i O'^SUM) 

IFCRESCAP LEO 0) GO TO 285 
IF(RESC(.P . GE TCA( IF) ) GO TO 286 
DO 831 J = i , IF 
DO 287 Jl=^i , IF 
IF < J1 . EQ . ) ) IFA = II- 
IF(Ji.GT 1) IFA-Jl-l 
IF(SCAP<IFA) FQ 0 0) L.O TO 287 
RESCAP-SGAP OFA)» ( 1 0-SUM) 

CNTEN = ACCM(IF,If A)/FLOAT(If ) 
ir<RESCAP or CN78N) GO 70 288 
ERnSN--=£RO*Ri:SCAP 
GO TO 289 

288 FROr>N--f pnJf.LNTFN 

289 ERS( IFA)^f RS( ir A)+Er.trif>N 
SUM=-SUM + EknSN/SCAP ( I F A ) 

IF(SUM CF 1.0) GO 70 2BG 

287 CONTINUE 
831 CONTINUL 
286 CONTINUfi 

DO 290 IFA=1,1F 

IF<SCAP<IfA) t(4 0 Oi r.(l TO 290 
EROSN-ERn*ACCM( IF , IFA 
ERS( IFA) =^F1<';( II A) - FROON 
SUM-SUM t f' ROSN/SCAP ( 1 f A ) 

290 CONIINUI 

284 CONTIMUt 

C UPDATE LOAD A1 NUDL PUINIS 

285 CONTINUE 
TFM-0 0 
liEDWN-0 0 

DO 2S0 If-=i ,Mf R 
CP Cl ( IC , ] ( ) 

WAT MAT ( IF; WAT MAT(1F ) *LRS(1f ) 
CCC=UAIMATi]f )/A( 

CAP=Sf:AH ( ]f )/AE 
XKP^CAP-0 y99»r:p 
XKOCAP 0 9V9»r.CC 
ADD-=0 . 0 

IF(CAP GT 0 0) AI)D-Cf I »CAP/ ( AF */KP*XKC) 

IF (ADD . L7 0 0 ) ADD-- 0 0 

FNLR = 1 0 fADD 

DELT--^DELX»F NLR/UF. L 

IF(DELT.LT DT5) GO TO 2400 

IFdC EQ INLS) GO TO 2S1 

CC]=CI(IC + 1 ,IF) 

CC( IC + 1 ,IF)=rCCI-*-(CCC-CCI )*DTS/I)FLT 
GO TO 2S2 

C SEDIMENT OUTfIOU FROM CHANNEL BLOCK 
251 IF(IT EQ.l) CO TO 253 

CCI-=G( IT-1 ,1 , IF)/a(IT-i , I) 

GO TO 254 

253 CCI=CI(IC,IF) 


.1.85 





254 G(IT,I,IK>=<CCI+<CCC-CCI)«DTS/DELT)*Q<IT,I) 

GO TO 252 

2400 XDEL=DELT*VEL/FNLR 
IFdC.EQ.i) GO TO 2401 
CCB=CC<IC,IF) 

GO TO 2402 

2401 CCB=CUP(IT,IF) 

2402 CC(IC + 1 ,IF)=CCEi+(CCC-CCBi)iKDELX/XDEL 

IF(IC EQ INLS) G<IT,I,IF)=CC<IC+1,IF)*Q(1T,1> 

252 TERS=-ERS( IF)*DTS/DELT 

IFCTERO.GT.BEDMAT(IF)) TERO=BEDHAT(IF) 

BEDMAT <IF)=BEDMAT(IF)-TERO 
TBM=TPM+EiEDMAT< IF) 

C UPDATE BED ELEUATION 

BEDWN=--f(EDWN+TERO/UDTH/POR (IF) 

250 CONTINIIF 

BEDFLU(IC)-BEDFLVCIC)-BEDUN 
C ADJUSThENT OF ACTIVE LAYER THICKNESS 
IF<TBM GT 0 0) GO TO 294 
DO 293 IF---1,NFR 

293 PCBC(IC;,ir)=PCBI(IC,IF) 

GO TO 204 

294 CONTINUE 
PARM^O 0 

DO 295 IF--=j ,NFR 

PCBC<IC,IF)=BEDMAT<IF)/TBM#100 0 
IF<SCAP<IF) GT 00) GO TO 295 
PARM = PARM't PCBCC IC , IF ) 

295 CONTINUE 

1F(PARH GI 0 0) GO TO 820 

DO 821 IFF=1 ,NFR 

IFA=NFR-IFF+1 

NARM=lf A 

DARM = DHM(If A) 

PARM^PCBC(IC.IFA) 

IF < PARM GT 0 0) GO TO 020 
821 CONTINUE 

820 CURARM-100 0/PARM*DARM/3Q4 8 
ADHT=^C;URARM-TBM/WDI M 
IF(ADHT IF 0 0) GO TO 204 
ADARM-ADHT 
TBM-0 0 

DO ?9fc IF - 1 .NFR 

BF.DMAT< IF ) i BEDMAT < IF ) t ADAKM*PCB I ( If, ,IF)/1 00 0*UDTH 

296 TBM = TBM<HE.DMAT( IF) 

DO 297 IF-1,NFR 

297 PCBCtIC,IF)^ BEDMAT <1F)/TBM*100 . 0 
GO TO 204 

213 CONTINUE 

C DEPOSITION CAUSED BY TF.RMINATION OF FLOW 



DO 260 IF=^1,NFR 
IFdC.EQ.INLS) GO TO 

261 

261 

CCdC+l,IF) = 0.0 

GO TO 262 

GdT, I, IF)= 0 . 0 


262 

BEDMAT (IF) ^'BEDMAT (IF 

)4CldC,IF)*AE 

260 

CONTINUE 


204 

CONTINUE 



RETURN 

END 
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SUBROUTINE DUBOYS<D,DFT,AE,WDTH,SLP,SPGR,SNU,GAMA,CHI, 

&UST,SFLOU) 

THIS SUBROUTINE COMPUTES THE TRANSPORT RATE OF NONCOHESIVE 
SEDIMENTS USING A DUBOYS TYPE FORMULA<GRAF, 1971) 

CALL SHIELDCDFT,SPCR,SNU,GAMA,UST,TC) 

CHI=CHI/<D)(<i|iO 75) 

TAO = GAMA*AE/UIDTH)*;SLP 
SFLOU=UDTH*:CHI*TAO»<TAO-TC) 

IF(SFLQW.LE.0.0)SrLOW=0 . 0 

RETURN 

END 

SUBROUTINE MEYER<MEYERP,USTAR,RUB,D,WE 1CHT,UTSC,S,UIDTH,Q,V,GB,XO) 

THIS SUBROUTINE COMPUTES THE TRANSPORT RATE OF NONCOHESIVE 
SEDIMENTS USING THE BEDLOAD FORMULA DEVELOPED BY MEYER-PETER 
AND MULLER(19^8) 

FCTN<F ,RO,REY)-< 1 /SORT < F ) )-2 , *AL OGl 0 (RD ) -1 . 1'1 + 2 . #ALOGi 0 ( 1 .+9 7S* 
iR0/(REY*5qRT<F) ) ) 

G=32.17 

REY-4 0*RHB*V/VI5C 

RO=2.*RHB/D 

RSTAR=2 . #l))«USTAR/VISC 

IF<RSTAR GT 70.0) GO TO 7 


F1=0 006 
F2=0 09 
1 = 0 

FCT1=FCTN<F1,RO,RFY) 

F3=0 S*(f l-'F2) 

FCT3=FCTNif 3,RO,REY< 

IF<FCTl»;t CT3) J ,5,2 
F2 = F3 
GO TO 3 
F1=F3 

1F(ABS(P2-F1) 0 0001)5,5,4 
I = I + 1 

IF < I L T 100) GO TO 6 
GO TO 7 
F-F3 
GO TO 0 

F=(l 0/(2 ♦AlOGjOCRinn 14>)*»2 
AK=V#SQRT(F/tt )/USTAF 
AK=0 75 

Y=(USTARiFUr,TAR )/( (S-1 0 )*G*D) 

G1=MEYERP*U1DTH*(UST AR**3.0)»(WEIGHT/G)»(S/(S - 1 .0 ) ) 
CALL SHIELD (D ,S , VISE. ,UE1GHT ,USTAR ,TC) 

C TSTAR = TC/( (S-1 . 0)H<WEIGHT*D) 

TSTAR=0.047 
G2=AK**1.5-(TSTAR/Y) 

IF(G2.LT.O.0) CO TO 9 
G2=G2**1.5 
GB=G1*G2 

XO=(GB»l,0E6)/(««WEICHT) 

GO TO 10 
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9 


10 

C 


06=0 0 
XO=0 0 
RETURN 
END 


SUBROUT 1NE YANG(DFT,UST,SNU,V,SLP,U,S,CONG) 


THIS SUBROUTINE COMPUTES THE TRANSPORT RATE OF NONCOHESIVE SEDIMENTS 
USING THE TOTAL LOAD FORMULA DEVELOPED BY YANG(1973) 


D=DFT 

A=UST»D/SNU 

IF(A GL 70.0) GO TO 7 

VCW=2 5/<ALOG10<A)-0.06)+0.66 

GO TO 8 

VCU=2.05 

ESP=U*SLP/U-VCW*SLP 
IF<E3P) 9,9,10 
CONC=0 0 
GO TO 11 

0 Fi=5.435-0 .2B6*ALOGi0(UiKD/SNU)-0 457*AL0G 1 0 ( UST/U ) 

F2=i.799-0 409»ALOG10 <W*D/SNU)-0 314»ALOGl0(UST/W) 

F3-ALOG10(ESP) 

E-F1+F2*F3 

c=io o»#r 

CONC=C 

1 CONTINUf- 
RETURN 
END 

SUBROUTINE SHIELD<0,S,VISC,WLICHT,USTAR,TC) 

THIS SUBROUTINE COMPUTES THE CRITICAL BED SHEAR STRESS DERIVED 
FROM SHIELDS' FUNCTION 
REY- UST AR*1)/VISC 
ir<RFY CT 10.0) CO TO 1 
TC--0 08*(S-1 0)*WCIGHT*D/REY*«0 4 
GO TO 3 

IF(RLY GT 500 0) GO TO 2 

TC'O n72*(S-l 0)»WEIGHT*D*RFY**0 16 

GO TO 3 

TC-0 n6*<S-l 0)*UFIGHT*D 

CONTINUE 

RE(URN 

END 

SUBROUTINE SETVFI.(D,UI) 

THIS SUHTTOLITINE COMPUTES THE SETTLING VELOCITY OF SEDIMENT PARTICLES 
DIMENSION A<2,11 ) 

DATA A<l,l),A<l,?),A<l,3),A(i,4),A(l,5),A(l,6),A(l,7),A(l,8), 
6A<1,9) ,A(1,10),A<1,11)/ 

60.04,006,0.10,020,040,080,150,200,3.00,700,10.00/ 

DATA A(2,i),A<2,2),A(2,3),A<2,4),A<2,S),A(2,6),A(2,7),A<2,B), 

4A < 2,9),A < 2,10),A(2,11 ) / 

60 14,0 32,0.76,2 20.5 30,10 50,16.90,20.30,25.60,39.50,44.00/ 
IF(D.LE.0,04) CD TO 20 
IFCD.GE.10.0) GO TO 21 
GO TO 22 

20 W=0.14/0.04#D 
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21 


GO TO 100 

W=4.5/3.0*(D-10.0)+44.0 
GO TO 100 
22 CONTINUE 

DO 10 1=1,11 

IF<ft<l,I) GT.D) GO TO 11 
Di=A<l,I) 

U1=A(2,I) 

GO TO 10 

11 D2=A(i,I) 

U2=A12,I) 

GO TO 12 

10 CONTINUE 

12 U=(U2-Ui)*<D-Di)/(D2-Dl)+Wl 

100 RETURN 

END 
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